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Abstract 


A real-time guidance scheme for the problem of maximizing the pay- 
load into orbit subject to the equations of motion for a rocket over a spheri- 
cal, nonrotating Earth is presented. An approximate optimal launch guidance 
law is developed based upon an asymptotic expansion of the Hamilton-Jacobi- 
Bellman or dynamic programming equation. The expansion is performed in 
terms of a small parameter, which is used to separate the dynamics of the 
problem into primary and perturbation dynamics. For the zeroth-order prob- 
lem the small parameter is set to zero and a closed-form solution to the zeroth- 
order expansion term of the Hamilton-Jacobi-Bellman equation is obtained. 
Higher-order terms of the expansion include the effects of the neglected pertur- 
bation dynamics. These higher-order terms are determined from the solution 
of first-order linear partial differential equations requiring only the evaluation 
of quadratures. This technique is preferred as a real-time on-line guidance 
scheme to alternative numerical iterative optimization schemes because of the 
unreliable convergence properties of these iterative guidance schemes and be- 
cause the quadratures needed for the approximate optimal guidance law can 
be performed rapidly and by parallel processing. Even if the approximate solu- 
tion is not nearly optimal, when using this technique the zeroth-order solution 
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always provides a path which satisfies the terminal constraints. Results for 
two-degree-of-freedom simulations are presented for the simplified problem of 
flight in the equatorial plane and compared to the guidance scheme generated 
by the shooting method which is an iterative second-order technique. 
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Chapter 1 


Introduction 


An approach to real-time optimal launch guidance is suggested here 
based upon an expansion of the Hamilton-Jacobi-Bellman or dynamic pro- 
gramming equation. In the past, singular perturbation theory has been used 
in expansion techniques used to solve optimization problems [1, 2, 3], For 
singular perturbation methods the states are split up into a set of ‘fast’ and 
‘slow’ variables. The solution is then sought in two separate regions; one re- 
gion where the fast states are dominant and an outer region where the slow 
states are determined. A composite solution can then be determined by com- 
bining the two solutions. Matching asymptotic expansions is one method for 
obtaining the final solution. This research uses a regular asymptotic expansion 
which is assumed valid over the entire trajectory of the launch optimization 
problem. An example of a launch optimal control problem is to determine the 
angle-of-attack profile which maximizes the payload into orbit subject to the 
dynamic constraints of a point mass model over a rotating spherical Earth. 
The solution of this type of optimization problem is obtained by an iterative 
optimization technique. Since the convergence rate of iterative techniques is 
difficult to quantify and convergence is difficult to prove, these schemes are not 
suggested to be used as the basis for an on-line real-time guidance law. 

In contrast, an approximation approach is developed which is based 
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upon the physics of the problem. Thrust and gravity are assumed to be the 
dominant forces encountered by the rocket while the angle-of-attack is usually 
kept small in order to minimize the effect of the aerodynamic forces acting 
on the vehicle. Numerical optimization studies [4] have been performed which 
support this assumption. These results also indicate that ignoring the aero- 
dynamic pitching moment has a negligible effect on the performance of the 
vehicle. Thus the launch problem would seem to lend itself to the use of per- 
turbation theory. It is shown that the forces in the equations of motion can be 
written as the sum of the dominant forces and the perturbation forces which 
are multiplied by a small parameter e, where e is the ratio of the atmospheric 
scale height to the radius of the Earth. The motivation for this decomposition 
is that for e = 0, the problem of maximizing the payload into orbit subject to 
the dynamics of a rocket in a vacuum over a flat Earth, is an integrable opti- 
mal control problem. The perturbation forcing terms in the dynamics produce 
a nonintegrable optimal control problem. However, since these perturbation 
forces enter in with a small parameter, an expansion technique is suggested 
based upon the Hamilton-Jacobi-Bellman equation. The expansion is made 
about the zeroth-order solution determined when e = 0. This zeroth-order 
problem is now solved routinely in the generalized guidance law for the Space 
Shuttle [5] with a predictor /corrector scheme employed to guide the vehicle 
along the desired path. 

The higher-order terms of the expansion are determined from the 
solution of first-order linear partial differential equations which require only 
integrations which are quadratures. Quadratures are integrals in which the in- 
tegrand is only a function of the independent variable. Previous solution meth- 
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ods applied to guidance problems have motivated the approach suggested here. 
These include the explicit guidance laws, E-guidance, developed by George 
Cherry [6] for the Apollo flight. By writing the dynamics strictly as functions 
of the independent variable a solution was obtained by quadrature integra- 
tions. Past applications [7, 8] of the proposed scheme, have shown that very 
close agreement with the numerical optimal path is obtained by including only 
the first-order term. Because no iterative technique is required, this scheme is 
suggested as a guidance law since the quadratures can be performed rapidly. 

Chapter 2 contains a general formulation of the perturbation prob- 
lem associated with the Hamilton-Jacobi-Bellman partial differential equation 
(HJB-PDE). The technique for determining the higher-order expansion terms 
due to the perturbation forces caused by the atmosphere and the spherical 
Earth model is discussed. Lastly, the recursive relationship for the control is 
presented. In Chapter 3, the characteristics for the Advanced Launch System 
(aka National Launch System) and the general equations of motion in terms of 
the small parameter e, are given. For e = 0, a simplified optimal launch problem 
in the equatorial plane is formulated, and its solution in terms of elementary 
functions is given in Chapter 4. The coordinate system transformation used 
to obtain the analytic solution is included. Also discussed is the linking of the 
trajectory subarc for the first stage to the subarc of the second stage. In Chap- 
ter 5 the first-order correction term to the control is determined. Results are 
presented in Chapter 6 and compared to the shooting method solution, which 
is a numerical iterative second-order optimization technique. It was found that 
during much of the first stage the aerodynamics are not small when flying the 
optimal vacuum trajectory. Chapter 7 presents a method for reshaping the 
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zeroth-order trajectory by including an aerodynamic effect. This effort cen- 
ters on the use of constant aerodynamic pulse functions which are obtained by 
averaging the aerodynamics along the zeroth-order path during various time 
intervals. Lastly, Chapter 8 relates perturbation theory and the Calculus of 
Variations with the expansion of the Hamilton-Jacobi-Bellman equation. The 
equivalence of the two solution methods is presented. 


Chapter 2 


The Peturbed Hamilton-Jacobi-Bellman Equation 


The optimal control problem can be formulated as one which mini- 
mizes a performance index subject to a set of nonlinear dynamics and a set of 
terminal constraints; that is, 

Minimize 

J = <t>(y/,T f ) (2.1) 

with the dynamics 

y = f(y,u,T)+eg(y,u,T) (2.2) 

subject to the terminal constraints 

^(2//> r /) = 0 (2.3) 

and the initial conditions 

y(t) = x = given (2.4) 

Note that y is an n-dimensional state vector, u is an m-dimensional control 
vector, e is a small parameter, r is the independent variable, y = dyjdr , t is 
the initial value of the independent variable, and x is the initial state at t. 

Eq. (2.2) is separated into two portions: primary and secondary dy- 
namics. Note that the control appears in both parts. The primary dynamics 
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can be assumed to dominate over the secondary dynamics because the sec- 
ondary dynamics are multiplied by the small parameter (e) and therefore have 
a small perturbing effect on the system. 

The Hamilton-Jacobi-Bellman (II-J-B) equation [9] is 

- P t = = min II = P x [r l + tg^ 1 } (2.5) 

u €.U 

where U is the class of piecewise continuous bounded controls and u opL {: r, P x , t ) 
is obtained from the optimality condition /7 U = 0 and from the assumption 
that the Legendre-Clebsch condition is satisfied (H uu is positive definite). In 
addition, f apt = f(x,u opt ,t ) and g ^ 1 = g{x ,u opt > 0- The Hamilton-Jacobi- 
Bellman equation will be used to determine the optimal control policy which 
minimizes the cost criterion J . 

The function P(x, t ) is called the optimal return function and is de- 
fined as the optimal value of the performance index for a path starting at x and 
t while satisfying the state equations (2.2) and the terminal constraints, i.e., 
P(x,t) = <t>(yf,T f ) at the hypersurface ^( y f ,r f ) = 0. The Hamilton-Jacobi- 
Bellman partial differentional equation (2.5) can be interpretated [10] as the 
derivative of the optimal return function P. The optimal return function is 
a constant since it is dependent only on the terminal conditions and thus the 
total derivative of the optimal return function along an extremal path must be 
zero. 

= Pt + Pr[/ op< + eg**} = 0 

at 

Each point in space belonging to the optimal trajectory must give the same 
value to the optimal return function as the optimal P{x,t) since the trajectory 
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is considered optimal from the initial conditions (x, t) to the terminal manifold. 
Now, if a non-optimal control is chosen at any point in the trajectory, then the 
resulting terminal state, as generated by the system equations, must produce a 
value for the optimal return function equal to or greater than the optimal value. 
Thus the control that minimizes the cost is the control which at each point of 
the trajectory causes the derivative of the optimal return function to be zero. 
This is the fundamental notion represented by the Hamilton- Jacobi- Bellman 
equation. Note that x and t can be either the initial or the current state and 
time, respectively. In this context, it will be used to represent the current state 
and time. Also note that every admissible trajectory must satisfy the terminal 
constraints = 0. 

P(x, t ) can be expanded as a series expansion in e as 

CC 

P(x,t) = J2 P dx,t)e' (2.6) 

1=0 

and the optimal control can also be expanded in a series expansion as 

OO 

u opt (x,P x ,t) = J2u i (x,t)t i (2.7) 

t=0 

where u opt is obtained by substituting Eq. (2.6) into Eq. (2.7) and expanding 
the function. Therefore, it is possible to obtain the control law in feedback 
form. 


The zeroth-order control , u Q) is the optimal control for the zeroth- 
order problem where c = 0. If an analytic solution can be obtained for the 
zeroth-order problem then higher-order solutions for the control can be ob- 
tained by expanding the Hamilton-Jacobi-Bellman equation 


/it 1 + r. 9it i 


t=0 


\i=0 


^i= 0 


i=l 


( 2 . 8 ) 
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where the dynamics have been expressed as expansions of the form 

OO 

r t (x,u" pt ,t) = J2Mx^t)e' (2.9) 

i=0 

oo 

g^ t (x,u opl i t) = , £9i(x,u,t)e i ( 2 . 10 ) 

t=0 

Expanding Eq. (2.8) and collecting terms of equal powers in e, produces the 
following set of linear, first-order, partial differential equations 

Pit + PijT = 1) 1 = 1,2,... (2.11) 

j= o 

= Pi_i,...,F 0 ) 

The expansion of the Hamilton-Jacobi-Bcllman equation will be detailed in the 
next section. 

2.1 Expansion of the H-J-B Equation 

The solution to the optimal control problem requires the evaluation of 
the Lagrange multiplier, P x . Note that the quantity P x is the partial derivative 
of the optimal return function with respect to the state y at the initial time 
or the current time (since at t = t, y = x). The function P x is expanded in 
a series in the small parameter e. The terms of this series expansion, P lx , are 
evaluated in terms of quadrature integrals which are functions of Ri . Recall that 
the functions R, require the previously evaluated terms P ]x , and 
for j = 1, ... ,i - 1. The coefficients /, and & are the i £/l term in the series 
expansion of / and g given in Eqs. (2.9)-(2.10). Since / and g are assumed to 
be sufficiently differentiable, they are expressible in a power series in e in terms 
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of the control. For a scalar control, this yields 

1 O'Ky, u, t) 


i=0 

oo 

g^(x,u^,t) = Y: 


t=0 


i! du ' 

1 d l 9{y,u,T) 


du 1 




( 2 . 12 ) 

(2.13) 


The above equations assume that the zeroth-order control, ito, is the dominant 
term in the series (Eq. (2.7)). This implies that the higher-order correction 
terms, u i, u 2 , ■ ■ •, have a much smaller effect on the optimal return function, 
P(x, <), than the zeroth-order term. The first four terms of / and g are obtained 
by use of Eqs. (2.12) and (2.13). 


/o = / <>pt (x,tio,0 = f(x,u 0 ,t) (2.14) 

/i = u 1 / u (x,a o ,0 (2.15) 

u 2 

y *2 = y/uu(x,Mo, 0 + U 2 f u (x,u 0 ,t) (2.16) 

u i 

h g /uuu(x, U.0, t) -f- U \ ll 2 fuu (x, ^0 j 0 

Tti3/ U (x, iia, t) (2.17) 


<7o = 

u 0 , t) = g(x,u Q ,t ) 

(2.18) 

9\ = 

u\g u {x,u 0 ,L) 

o 

(2.19) 

92 - 

Ul 

yf7u u(x, u 0 , l) + u 2 g u (x, Uq, t ) 

(2.20) 

93 = 

u\ 

~4 9 uuu(x , n 0 , 0 + Uiu 2 g uu (x, u 0 , t) 
o 



-Hx 3 £ u (x, it 0 ,i) 

(2.21) 
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Note that in taking the partials with respect to u in Eqs. (2.12) and (2.13), the 
partial is taken first and then the partial is evaluated at x, t with e set equal to 
zero. In other words, the partials are evaluated along the zeroth-order path. 


2.2 Solution by the Method of Characteristics 

The H-J-B equation (Eq. (2.5)) is a first-order partial differential 
equation. The expansion of the H-J-B equation results in the first-order dif- 
ferential equation for P, stated in Eq. (2.11) with the boundary condition 
Pi(xf,t f ) = 0, for i = l,---. Recall that j^ pt denotes the dynamics of the 
zeroth-order problem (e = 0) using the zeroth-order control u = Uq. Recall also 
that the forcing term P* is only a function of expansion terms of P of order 
less than i. 

The method of characteristics is used to solve a set of linear or quasi- 
linear partial differential equations. This technique [11] requires the identifi- 
cation and solution of characteristics curves. The characteristic direction ds is 
defined by the equation 

P it (dr), 4 - P ix ( dy), = (dPi), *= 1 , 2 ,--- ( 2 . 22 ) 


Eqs. (2.11) along with (2.22) can be put in the form 


1 /o 


Ri 

(dPi) a 


(2.23) 


10 dr). (dy)s\[Pu\ L {dPt) s \ 

The characteristic directions for Eq. (2.23) are given by the solution of the 
differential equation that is obtained by setting the determinant of the matrix 
given in Eq. (2.23) equal to zero, such that 


(dy), - fo(dr), = 0 


(dy/dr) s = f 0 


(2.24) 
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The subscript s denotes the characteristic direction. Therefore, the charac- 
teristic curves of the equations, for any order term of P iy are given by the 
zeroth-order optimal trajectory 

Vo — fo (2.25) 

whose solution is denoted as y 0 (r;x, t). 

The solution for P, is given by 

Pi(x,t) = - J h n°dr (2.26) 

where is defined along the zeroth-order path as 

Ri = ^( 2 /o,r, Pi-i(yo,T), -,P 0 {y 0) T)), * = 1,2,--- (2.27) 

Therefore, having already determined P terms of order less than i, a solution 
for Pi can be determined by integrating Ri from the current ‘time’ to the final 
‘time’ along the zeroth-order path. 

2.3 Determination of the Optimal Control 

Since the primary and secondary dynamics, / and g , are expanded 
in terms of the control (Eqs. (2.12) and (2.13)), the control expansion terms 
ito, Ui, v. 2 , • • •, need to be determined. The optimality condition provides the 
necessary tool to obtain these control terms. It can be stated as 

«,[/« + €9.)= E p .. c ' [£(/.. + £ 9.-K =0 (2.28) 

L i=0 J t=0 -I 

By expanding and multiplying out the terms of the two power series and equat- 
ing like powers of e, the following relations are obtained 


P Ox fu = 0 


(2.29) 
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e 1 

P ) r [,9u + *i|/uu] + P\ x fu = 0 

(2.30) 

e 2 

• Po z 1 9uu 4“ -Uifuuu 4" ^ 2 /uuj 



+ P\ z {gu + U\fuu] + Plxfu = 0 

(2.31) 


Note that Uo, the optimal control for the zeroth-order problem, can 
be solved using Eq. (2.29). Similarly, u x can be solved using Eq. (2.30) and u 2 
can be solved using Eq. (2.31). 

2.4 Determination of the Forcing Functions 

Eqs. (2.14) — (2.21) and (2.29) — (2.31) can be used to solve for the 
forcing functions ft, where Eq. (2.11) can be restated as 

ft« = — Pjz(fi-j 4~ 9i- j-i) z=l,2,... (2.32) 

j = o 

Using the above equations, fti is 

fti = —Po x (f\ + 9o) = —Po I (uifu 4- <j) (2.33) 

With the use of the optimality condition of Eq. (2.29), fti becomes 

ft. = -Po x 9o (2-34) 

Similarly, the equation for ft 2 is 

ft 2 = — Po x (f 2 + fti ) — P\ x (fi + 9o) (2.35) 

ft 2 simplifies to the following equation when Eqs. (2.14) — (2.21) and (2.29) — (2.30) 
are substituted into the previous equation. 

ft 2 = yftu/u U - P. x flb 


(2.36) 
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Finally, R 3 can be expressed as 


R-3 = —Po z (h + . 72 ) — P \ x (h + <7i) — P 2 x(/i + go) 


(2.37) 


This simplifies to 


Ra = 


Po x 
= Pox 


* 4 , 

3 Su 


u 


12 


If 

rtJ U 


+ U\U2fuu ~ 



U\U 2 . 
T o 


■Pi. 


2 

uu\+Pxx[ U ±U\-P2 I 9o 

[-^■<7u] - p2x [50 + ^/u] (2.38) 


Using the expression for P, the expression for the Lagrange multipli- 
ers, P tz , can be expressed as 


Pix = 


dp 

dx 


r l f dp , nl ut „ . uif 

= -/ + 


dt 1 


(2.39) 


Once these Pi z are determined, they can be used in the optimal control ex- 
pansion (Eq. (2.7)). As made apparent in the above equations, the solution 
becomes increasingly complex as the higher-order correction terms rely on the 
state information from the lower-order trajectories. 



Chapter 3 


Modelling of the ALS Configuration 


This chapter presents the modelling characteristics and the equations 
of motion for the rocket. Included are sections on the properties of the propul- 
sion, aerodynamics, masses, gravity, and the atmosphere. A small expansion 
parameter, the ratio of the atmospheric scale height to the radius of the Earth, 
is then used to separate the dynamics into the primary and perturbation ef- 
fects. Lastly, the equations of motion for the zeroth-order problem of flight in 
a vacuum over a flat Earth are presented. 

The Advanced Launch System (ALS) is designed to be an all-weather, 
unmanned, two-stage launch vehicle for placing medium payloads into a low 
Earth orbit. The spacecraft (fig. 3.1) consists of a liquid rocket booster with 
seven engines and a core vehicle that contains three engines. All ten liquid 
hydrogen/liquid oxygen low cost engines arc ignited at launch. Staging occurs 
when the booster’s seven engines have exhausted their propellant. The three 
core engines burn continuously from launch until they are shut down at or- 
bital insertion. Launched in the equatorial plane and ending at the perigee 
of a 80nm by 150nm transfer orbit, the flight occurs in two-dimensions over a 
nonrotating, spherical Earth. Note, the booster is assumed to ride on top of 
the core throughout the first stage trajectory. 


14 



16 


3.1 Equations of Motion for the Launch Problem 

The general equations of motion for a launch vehicle modelled as a 
point mass over a spherical, nonrotating Earth are given for flight in three- 
dimensions as 


h = 
V = 

V sin 7 

(T cos a cos 0 — D) 

v — n cm "Y 

(3.1) 

(3.2) 

— — y oui j 

m 

[-(T cos asin0 - Q) sin /x + (T sin a + L) cosp] 

7 = 

mV 



r V 9 i 

+ — cos 7 

+[ (r e +h) E J 

(3.3) 


\{T cos a sin 0 — Q) cos fx + (T sin a 4- L) sin /r] 


X = 

(mV cos 7) 



V tan 0 cos 7 cos x 
+ (r e + h) 

(3.4) 

6 = 

V cos 7 cos x 
(r e + h) cos 0 

(3.5) 

4> = 

V cos 7 sin x 
(r e -1- h) 

(3.6) 

m = 

o' T vac 

(3.7) 


The vehicle coordinate system is shown in figure 3.2. Note, the engines are not 
gimbaled and the aerodynamic pitching moments are neglected. For a vertical 
launch Eqs. (3.3)-(3.4) experience a singularity caused by the velocity being 
zero and by a flight path angle of 90 degrees, respectively. Therefore, a pitch- 
over maneuver must be made at launch and equations of motion written in a 
different coordinate frame must be used. 
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Figure 3.2: Coordinate Axis Definition 
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3.2 Propulsion 

Thrust is assumed to act along the centerline of the booster-core vehi- 
cle configuration and to be the same constant value for each engine. The total 
thrust of the rocket changes after staging as the seven engines of the booster 
are discarded, leaving only the three engines of the core vehicle. 


T = (T vac - npA e ) T vac = n x 580, 1 10. lbs. 


where T vac is the total value of the thrust when acting in a vacuum and the 
number of engines is n — 10 for the first stage and n = 3 for the second stage. 
Notice the variation of the thrust due to the atmospheric pressure p is given 
for an underexpanded nozzle and thus a conservative value for thrust is used. 
The value of the engine nozzle exit area is A e = 5814.8/144. sq ft. The specific 
fuel consumption of the rocket is 


and the specific impulse I sp 
after staging occurs. 


sec 


I sp i] S ft 


(3.8) 


= 430. seconds. The value of a remains the same 


3.3 Aerodynamics 

Since sideslip causes drag, the vehicle is assumed to fly at zero sideslip 
angle, so that only the angle-of-attack gives the orientation of the vehicle rel- 
ative to the free stream. The direction of the lift vector is then controlled 
through the velocity roll angle. With no sideslip, the side force Q is identically 
zero. Therefore, 
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L = C L qS , D = CnqS, Q = C Q(? S = 0 (3.9) 

where Cx, Cd,Cq are the lift, drag, and side force coefficients, respectively, S is 
the cross-sectional area of the combined vehicle (booster 4- core), and q = \pV 2 
is the dynamic pressure. The cross-sectional area S is assumed to be the same 
constant value before and after staging occurs. 

The aerodynamic data has been provided in tabular form [4] and is 
modelled by polynomials in a with Mach-number-dependent coefficients. For 
the first stage, the aerodynamic coefficients are written as 

C D (M,a) = C Do (M) + Co a2 (M)a 2 + C DQ3 (M)a 3 
C L (M } a) = C La (M)a (3.10) 

where the Mach-number-dependent terms have been obtained from cubic-spline 
curve fits of the tabular data. Three-dimensional plots [12] of the first stage 
drag and lift models are shown in Figures 3.3 and 3.4. Note that the drag 
coefficient of this vehicle at supersonic and hypersonic speeds has a minimum 
at a positive angle of attack as shown in Figure 3.3. This is caused by the 
aerodynamic shielding of the booster by the flow field of the core. 

After staging, the vehicle operates in the hypersonic flow regime and 
the aerodynamic force coefficients are modelled as 

Cd(cc) = Co a + Co a & + Co a 2 & 2 

C L {a) = C La a + C La y (3.11) 

with constant coefficients Cq 0 = .2011, Co a = 0.0, Co a2 = 001811, Ci a = 


21 



Figure 3.5: Second Stage Aerodynamic Model 

.039962, and C^ q2 = .00100272. The aerodynamic plot of Ci and Co is pro- 
vided in figure 3.5. 


3.4 Mass Characteristics 


Ihe inert weights of the booster and core, the weight of the propellant, 
the payload and payload margin, and the weight of the payload fairing comprise 
the ALS takeoff weight. The fairing encases the payload and is carried along by 
the core vehicle until orbital insertion. The vehicle mass and sea-level weight 
characteristics are shown in Table 3.1. The time at which staging is to occur is 
obtained from the first stage mass flow rate and the propellant of the booster 


t 


stage — 


m propellant 

7aT vac 


153.54 sec. 
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Vehicle Stage 

Vehicle Component 

Take-off Weight 
(lbs.) 

Core 

Inert Mass 

176,130.00 

Propellant 

1,479,180.00 

Payload 

120,000.00 

Payload Margin 

' 12,000.00 

Payload Faring 

39, 120.00 

Total Core 

1,826,430.00 

Booster 

Inert Mass 

216,880.00 

Propellant 

1,449,980.00 

Total Booster 

1,666,860.00 

Core + Booster 

Total at Take-ofT 

3,493,290.00 


Table 3.1: Vehicle Mass Characteristics 


where the vacuum thrust per engine is Tya c 580110* 

Once the stage time, the total first stage mass flow rate, the takeoff 
weight, and the inert weight of the booster are known, then the weight of the 
vehicle at the end of the first stage and the initial weight in the second stage 
can be calculated. For this vehicle the values are 

m S iage\ = 1421890. lbs., m stagc2 = 1250010. lbs., A m stage = 216880. lbs. 


3.5 Gravitational and Atmospheric Models 

The gravitational acceleration is modelled as an altitude-varying func- 
tion by the inverse square law, 


23 


but will be assumed constant in the zeroth-ordcr problem to facilitate obtaining 
an analytic solution. The constant values for gravity at sea-level and for the 
radius of the Earth are 

g s = 32.174 — - r e = 2.09256725 x 10 7 ft. 

sec** 

The atmospheric density is expressed by the exponential function, 

p = p r e~ {re+h)/h ‘ = p r e~ T * /h ‘e~ h/h ’ - p s e~ h/h ’ (3.12) 


where h a is the atmospheric scale height and p s is the sea-level reference density. 
The values for these parameters are 

p s = .002377 h, = 23, 800. ft. 

ft J 


The form of the density is chosen to motivate the selection of a small 
parameter to exclude the aerodynamics in the zeroth-order dynamics. If e is 
chosen as 

e = h s /r e (3.13) 


and defining 

6(e, h) = (3.14) 

t 

then by atmospheric properties <5(e, h) > 0. The exponential density also sat- 
isfies the requirement [3] that the perturbation term in the dynamics remains 
small, i.e., 

lim <5(e, h) — + 0 (3.15) 

Satisfaction of this property will allow more general atmospheric models to be 
used in the launch problem. 
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The atmospheric pressure is also expressed as an exponential function, 

p = p 3 e~ h/>Xp (3.16) 

where hp is the atmospheric pressure scale height and p 3 is the sea-level reference 
pressure. The values for these parameters are 

lbs 

p s = 2116.24^ hy = 23,200. ft. 

The speed of sound can be obtained by the relationship 

/TP 

SOS = A — 

V p 

with the specific heat ratio for air given as T = 1.4 . 

The gravity can be rewritten as 


g 3 h( 2r e -f- h) eg s h(2r e T h)r e 

9 — 9s 7 . i \ o — 9» ~ 


(3.17) 


(r e + h) 2 JS h 3 (r e + h ) 2 
where the expansion parameter has formally been introduced and the second 
term is clearly small in comparison to the first term which is the value for 
gravity at sea-level, g 3 . 


3.6 Expansion Dynamics 

In terms of the small parameter e, the full-order equations of motion 
are rewritten as 


h 

V 


V sin 7 
71 


m 

+e | 


cos a cos j3 — g 3 sin 7 

npA e r e „ . g 3 h(2r e + h)r e sin 7 

cos a cos p + 


mh 3 


h 3 (r e + h) 2 


(3.18) 

(3.19) 

pSV 2 C D r c ' 

2mh s 
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T 

*■ in 




— e 


mV 
npA e r e 

mV h. 


(cos a sin /3 sin fi - sin a. cos p) 


fh cos 7 


V 


'cos a sin j3 sin p — sin a cos p) 


pSVr e . _ . 

+e— — ; — (C q sin fi + G/, cos fi) 


+e 


2mh s 

V 


+ 9. 


h(2r e + h ) \ r e 


r e + h V ( r c + h) 2 I h. 


cos 7 


(3.20) 


T 


mV cos 7 
np/ l e r, 


+e 


mVh s cos 7 
Tnh s cos 7 


(cos a sin cos fi + sin a sin p) 

(cos a sin (5 cos fi + sin a sin fi) 


{Ci, sin p — Cq cos p) 4- 


V r,, tan O cos 7 cos \ 
h s (r e + h) 


■ V cos 7 cos v h . 

0 = 

r e cos (p a „ 

• V cos 7 sin * /i 

0 = ^ — (l - e h? 


(3.21) 

(3.22) 

(3.23) 


Where the binomial formula has been used to rewrite (r e +k) 1 for the longitude 
and latitude since r e h. 


3.6.1 Two-Dimensional Flight 


In this section the three-dimensional equations of motion are reduced 
for flight in a great-circle plane (the X-Z plane) over a flat, nonrotating Earth. 
If the vehicle is assumed to be restricted to fly in the equatorial plane then 
the lift, thrust, and velocity vectors all lie in the same plane and the roll angle 
(p = 0) is eliminated from the equations. Under the previously mentioned 
assumptions of no side force ( Q = 0 ) and no sideslip (/? = 0 ), the zeroth-order 
equations of motion representing flight in a vacuum over a flat Earth become 


h = 


(3.24) 


V sin 7 
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V 

7 

9 

m 


cos or — q s sin 7 

m 

Tvac . Qs 

— — sin a — — cos 7 
mV V 

V cos 7 




m = m 0 - aT vac {r - r 0 ) 


X = Xo = 0.0 


<(> = = 0.0 


(3.25) 

(3.26) 

(3.27) 

(3.28) 


These are the system dynamics used to obtain an analytic solution to the 
zeroth-order optimization problem presented in the next chapter. 


Chapter 4 


Zeroth-Order Optimization Problem 


The solution to the zeroth-order optimization problem is derived by 
a coordinate transformation. A canonical transformation from the wind axis 
to the rectangular or local horizon coordinate frame allows the zeroth-order 
problem to be solved analytically. The solution is in closed form up to some 
constants that can be determined numerically to solve the two-point boundary 
value problem. The conditions for connecting the second stage subarc to the 
first stage subarc are then presented. 

4.1 Optimization Problem Statement 

In this section the zeroth-order optimization problem is presented. 
The problem is to maximize the payload into orbit 

J = —m/ 

subject to terminal constraints on the altitude, velocity, and flight path angle, 
h f = h U P ^ V f = v /s pec > 7/ = 7/ tpcc 

subject to the state discontinuity in the mass at a interior point where staging 
occurs, 

Wastage 2 stage 1 stage 
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and subject to the equations of motion for flight in the equatorial plane. 


h = 

V sin 7 

( 4 . 1 ) 

V = 

T 

— cos a — g s sin 7 
m 

( 4 . 2 ) 

7 = 

T . g s 

— — sma- 77 cos 7 

mV V 

( 4 . 3 ) 

0 = 

V cos 7 

r e 

( 4 . 4 ) 

m = 

—oT =>• m = m 0 — oT(t — r 0 ) 

( 4 . 5 ) 


Note, in this section and when discussing the zeroth-order trajectory, the total 
vacuum thrust will be represented by T and the subscript notation will be 
dropped. 

The Hamiltonian for this system can then be expressed as 

T T g s 

H = X h V sin 7 + X v { — cos a - g s sin 7) + A^( — - sin a - ~ cos 7) (4.6) 

m mV v 

The zeroth-order control law determined by the optimality conditon is 

T T 

H a = AvsinaH — A~cosa; = 0 (4.7) 

m mV 


By the strengthened Legendre-Clebesch condition H aa > 0 choose 


tan a 
cos a 

sin a 


A_ 

VX v 

V Xy 

“^A v ) 2 + y 

^2 

+ y 


( 4 . 8 ) 


Whereas the optimal control can be derived in terms of the states and Lagrange 
multipliers, an analytic solution is not possible for the states and Lagrange 
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multipliers written in the wind axis frame. Therefore, a coordinate transfor- 
mation into the Cartesian reference frame is presented in the next section. In 
section 4.3 an analytic solution is obtained using this transformation. 

4.2 Zeroth-Order Coordinate Transformation 

The analytic solution for the zeroth-order problem can be found in 
the Cartesian coordinate system but the equations of motion of the full sys- 
tem which include the aerodynamic forces are written in the wind axis system. 
Therefore, to derive the zeroth-order control and the first-order correction to 
the control the transformation of coordinates and especially the transformation 
of the Lagrange multipliers must be known. This can be accomplished by a 
canonical transformation [see appendix B] from the (0, <f > , h) coordinates to the 
right-handed coordinate system (X, Y, Z), where X is positive in an eastward 
direction along the equator, Z is positive pointing towards the Earth, and Y 
is orthogonal to the X — Z plane. The relationship between the two reference 
frames (see figure 4.2) is X = r e 6, Y = r e 4 > , and Z = —h. In two-dimensions, 
the corresponding velocity coordinates (u,w) are considered positive in the pos- 
itive X and Z directions, respectively. A necessary and sufficient condition [13] 
for a canonical transformation is the equivalence of the Hamiltonians in the two 
reference frames. 


Hlh = A xdX -|- XydY + A hdh + A U du + A W dw 
Hwind — XodQ + \<t,d(j> + A hdh + A vdV -I- X-jd'y 


(4.9) 

(4.10) 



T 



T W 



Figure 4.1: Transformation of Coordinate Systems 


31 


This equivalence is obtained through the Jacobian of the transformation. There- 
fore, the transformation 


u = V cos 7, w = — V sin 7 


requires 


Av 


"1 

du 

dw 


' A u ' 

Ajy 

= 

& 

ft 


J 

. ^7 




and thus, 




cos 7 

— sin 7 

' A u ‘ 

. . 


— Esin 7 

— V cos 7 

Ayj 


This produces the transformation of the Lagrange multipliers, 


Av = A u cos 7 — Au, sin 7 
A 7 = -l / (A u sin7 + AujCOS7) 
A 0 = f e Xx 

X4, r e Ay' 


and the transformation of the states, 


V = Vu 2 + ui 2 
w 

sin 7 = 

1 V 


(4.11) 


(4.12) 

(4.13) 

(4.14) 

(4.15) 


(4.16) 

(4.17) 


4.3 Zeroth-Order Analytic Solution in the Cartesian 
Frame 

In this section an analytic solution will be derived for the zeroth-order 
problem of maximum payload into orbit for flight in a vacuum over a flat Earth. 
This solution is made possible by the coordinate transformation presented in 



32 


the previous section. The equations of motion in a Cartesian coordinate frame 
are 


X = 

n 

(4.18) 

Y = 

0 
II 

II 

1 
o 


h = 

—w 

(4.19) 

u = 

T 

— cos Op 
m 

(4.20) 

V = 

0 
II 

§ 

II 

1 
o 


w = 

T 

sin Op + g 3 

m 

(4.21) 

m = 

—oT ==» m = m 0 — <jT(t — r 0 ) 

(4.22) 




( 4 . 25 ) 
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The Lagrange multipliers are obtained using X y = -//J 


Ax = 0 

A/i = 0 

A u = —Ax 

X w — Xh, 


with the boundary conditions 


Ax(t/) = 1'x, A h {Tf) = U h) 


X u (Tf) = IS U , X w {tj) = v x 


where i 'x,Vh,v u , v w are unknown Lagrange multipliers associated with the ter- 
minal constraints. For the unconstrained downrange problem, the solutions to 
the adjoint differential equations are 


Vx = o 


Vh 

(4.26) 

t'u Cy 

(4.27) 

c w + X h (r - To) 

(4.28) 


The equations of motion can be integrated by changing the independent vari- 
able from time to mass and using the mass equation (Eq. (4.5)) to substitute 
mass for r. As a consequence, the Lagrange multipliers are rewritten as 


A u 
A*, j 


A’ + A’ 


= C u 

_ 7* _ \ — 

(jT 

= cm? + bm + a 


(4.29) 

(4.30) 

(4.31) 
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where 


A 2 , 

c wry 

(4.32) 

b — — ^ 

ai 

(4.33) 

a = Cl+Vl 

(4.34) 

C w = C„ + 

(4.35) 


The derivatives of the states with respect to mass are 


du 

dm 

dw 

dm 

dX_ 

dm 

dh_ 

dm 


am/ cm 2 + bm + a 
A w 


am/crn 2 + bm + a aT 

u 

~Vr 

w 

aT 


(4.36) 

(4.37) 

(4.38) 

(4.39) 


Note that c > 0, a > 0, and the discriminant of the quadratic mass equation 
A = 4ac — 6 2 > 0 since 


A = 


(ary 


(a h c u y 


(4.40) 


From these differential equations the solution is found from standard integrals. 


U — Uq — 


W — Wq 


c . u 

a /a [ 
9s_ 

‘ aT 

\h 


sinh 1 


'2a + bm 

7T 


m 


^ — sinh 1 ^ 


_ i ( 2a + brrio ' 
mo-v/A , 


(4.41) 


a 2 T y/c [ 

c, 


(m - m 0 ) 

sinh- 1 ( 2cm Z - ) - sinh- 1 


a /a [' 

h = h o — 


sinh 1 
9s 


( 2 cm + b\ 

\7E~) 

2a + 6m )-sinh-‘ ( 


y/E 


V VA 

_ j ( 2 a + bmo \ 
mo'/K /. 


(4.42) 


2 (aTy 


, \2 , (m-mo) 

(m - m 0 ) + — — — 


Wq 
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f, 9s , \2 , %) 

*» " 2My' m - m °> + - aT w ° 


" 2(crT) 2 u ' 

i ] 

~ m 7 7= sinh 


-ra - - . - sinh 

cr((jy jv/a 


h _,/2cm+_6 

V v'S 

_, (2a + bm\ 


_ I / 2c7Tio 4- b 

. i ( 2 a T buiQ ^ 


cr(aT) 2 c 


'em* 4- bm 0 + a - V cm? + bm + 


At the final time, /// = - 1 by the transvcrsality condition. Using the Hamil- 
tonian and the three state equations u,w, and h, which have prescribed initial 
and final values, the four unknown constants associated with the two-point 
boundary value problem can be solved. For the problem of flight restricted to 
a plane, the unknowns are mf,C u ,C w , and A/,. The analytic state equations 
(Eq. (4.41)-(4.43)) are nonlinear and thus no statement can be made about the 
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existence or uniqueness of the set of constants found. Therefore, if multiple 
solutions are found the solution set which minimizes the Hamiltonian would 
be chosen. At the very least, the Legendre- Clebesch condition, H uu > 0, for a 
weak relative minimum must be satisfied. 


4.4 Linking the First and Second Stage Subarcs 

Of interest in this section is the linking of the two subarcs of the 
two-stage rocket. By the corner conditions, the Lagrange multipliers for all the 
states must be continuous. 

Ay(Ls(age — ) = A y (t s iage + ) (4.44) 

The analytic solution previously presented is still valid for either subarc but 
only by using this relationship between the Lagrange multipliers can the sec- 
ond stage be connected to the first stage subarc. Recall that the constant C w 
is associated with the initial condition of the Lagrange multiplier for the ver- 
tical velocity component. For a subarc with first stage initial conditions, the 
equations become 

^w(tstage) = Gui + ^h(tstage ~ to) (4.45) 

A u ,(f) — A wastage) + A /,(£ — tstage) t > fstageT (4.46) 

Rewriting the Lagrange multipliers using the corner condition and with mass 
replacing time as the independent variable, results in 

\ h = v h — constant to < t <tj 

A u = 


u u = C u = constant to < t <tf 


(4.47) 

(4.48) 
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\v = c w + -m) to < t < t stage 


(4.49) 


Atu(M-) = C w + -£r(m 0 - m st a Q c\) 4- -~(m stage2 - m) (4.50) 

(71 l (71 2 

where T\ and T 2 represent the thrust for the first and second stages, respectively. 
The equations of motion, written with mass as the independent variable, which 
were previously presented are still valid but the constant coefficients of the 


quadratic equation are of a different form. 

A„ 4- A^ = c'm 2 + b'm + a 


b' = -4-^. 

(71 2 

a = Cl + C'l 


stage 2 (dstagc\ 


c' w = c w c A,^ + A,(^F- 

(71 i (7 L 2 

7 ^ . x / stage 2 TI7stage\ \ 

- + — — — ) 

(712 (71 1 


Therefore, the state equations become 


Cu 

amydrrd 4- b'm + a' 

Au> Qs 

amy/dm 2 4- b'm 4- a' (tT 2 


(4.51) 

(4.52) 

(4.53) 

(4.54) 

(4.55) 


The same standard integrals apply to the solution of the problem because 
a 1 > 0, d > 0 and the discriminant 


A' = 4aV — 5' 2 = 4 ( Cu>0- 


(4.56) 
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The simplified form of the solution to the state equations (Eqs. (4.41)-(4.43)) 
is also still valid but with the first stage subarc used as the initial conditions 
of the second stage subarc. 
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These are the equations that result from linking the first stage subarc to the 
second stage subarc. These equations will be used to evaluate the states at a 
time after staging occurs when the initial time is before staging. The first-order 
correction terms will require the analytic solution for the states at any future 
time along the zeroth-order trajectory. 



Chapter 5 


First-Order Corrections 


The use of the asymptotic expansion of the dynamic programming 
equation as discussed in Chapter 2 by the approximate optimal guidance scheme 
is an improvement over past analytic techniques whose guidance laws were lim- 
ited to operate in the exoatmospheric region [6, 14]. The higher-order correc- 
tion terms of the HJB expansion can be used to compensate for the effects of 
the atmospheric forces neglected in the exoatmospheric solution. The deter- 
mination of the first-order correction to the zeroth-order control is the subject 
of this chapter. As noted before, the solution to the first-order optimization 
problem requires only the integration of quadratures, which can be evaluated 
quickly enough to permit this method to be implemented as a real-time guid- 
ance scheme. The correction to the Lagrange multipliers and thus the cor- 
rection to the control is constructed in the following sections. Also derived 
are all the partial derivatives needed to evaluate the quadratures. The partial 
derivative chain rule is employed since the analytic solution is found in the 
Cartesian frame while the first-order forcing function, R\, used to evaluate the 
quadratures is expressed in the wind axis frame. Recall that the angle-of-attack 
is the control variable and the aerodynamic coefficients are modelled as func- 
tions of the angle-of-attack. For this reason the perturbation dynamics are left 
expressed in the wind axes frame. 


40 


41 


5.1 Correction to the Lagrange Multipliers 


The higher-order terms of the optimal return function were presented 
in Eq. (2.26). 


= -J t/ n°dr 

By taking the partial derivative of this integral the correction term to the 
Lagrange multiplier can be calculated. Recall, 
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where the first-order forcing function was R\ = —Po x <)a- 

The first-order correction term for the Lagrange multipliers is used 
to determine the first-order expansion term of the control. By the first-order 
optimality condition, Eq. (2.30), the correction to the control is obtained. 


u \ — — (fuuPo x ) ' [Po z 9u T P u/u. 
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5.2 The First-Order Forcing Function 

For the launch problem as formulated in the wind axis frame, the 
first-order forcing function is 
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The Lagrange multiplier for the first-order term of the expansion series is found 
by integrating the partial derivative of R\ with respect to the initial state. For 
the launch problem, the optimal control depends on the Lagrange multipliers 
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for the velocity and flight path angle, i.e., x = [Vo, 7 o]- The partial derivative 
of the first-order forcing function with respect to the initial state is 
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The paxtiaJs of the wind axis states and Lagrange multipliers are related to 
the partials of the analytic Cartesian states and Lagrange multipliers by the 
canonical coordinate transformation. These partial derivatives are presented 
in subsequent sections. 


5.3 Relating the Partial Derivatives of the Wind Axis 
Frame to the Partial Derivatives of the Cartesian 
Frame 


The canonical transformation of section 4.2 provides all the infor- 
mation needed to relate the analytic solution of the zeroth-order states and 
Lagrange multipliers to the states and Lagrange multipliers in the wind axis 
frame. Thus, the variations in the analytic Cartesian coordinates due to varia- 
tions in the initial wind axis states can be determined and it was for this very 
reason the canonical transformation was necessary. Using the relationships 
obtained in section 4.2, the partial derivatives of the wind axis coordinates 
become 
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and from the zeroth-order control law Fq. (4.8) 
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Now that the partial derivatives for the wind axis coordinates are expressed 
in terms of the partial derivatives of the Cartesian coordinates, the partial 
derivatives of the Cartesian coordinates with respect to the initial states are to 
be derived along the analytic zeroth-order trajectory. 


5.4 Partial Derivatives of the Analytic Solution 

In this section, the partial derivatives of the Cartesian coordinates are 
derived. The zeroth-order analytic trajectory is used to evaluate the integral of 
the partial of the forcing function R] from the initial time to the final time. For 
the sake of notational brevity, the following common terms and their partial 
derivatives are defined. 


5.4.1 Partial Derivatives of Some Common Terms 

The partial derivatives of the constants a,b,c, and C w used to express 
the analytic state equations arc 
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dc 

dx 
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Recall that the function A = 4ac — b 2 , so the partial derivative is 
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Let the arguments of the inverse hyperbolic sine function be denoted 
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and by the partial derivative chain nile for a trignometric function, the partials 
of the inverse hyperbolic sine functions are 
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5-4.2 Partial Derivatives of the Analytic States 

The general form of the state equations in Eqs. (4.41)-(4.43) is used 
to derive the partial derivatives of the states with respect to the initial veloc- 
ity or flight path angle. Using the terms defined in the previous section and 
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simplifying the equations, the partial derivatives arc 
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The initial velocity components expressed in terms of the wind axis states are 


Uo = V 0 cos 70 , w 0 = - 1 / 0 sin 7 q 


(5.27) 


and therefore the partial derivatives with respect to the initial velocity and 
flight path angle are 
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These partial derivatives are valid for a point during the first or second 
stage of the trajectory with initial condition corresponding to that subarc. For a 
point on the second subarc with first stage initial conditions, the state equations 
which link the two subarcs must be used. Note also that these equations all 
depend on the partial derivatives of the constants, A h ,C u ,C w , and m f which 
are unknown. The partial derivatives of the constants are dependent on the 
initial and final conditions of the two-point boundary value problem. Using 
the transversality condition 
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r r /dhf &Uf dwf 9 Hr\ i. 

These results produce a system of four equations (-g? , "a* . “5*V lineax 


in the four unknown partial derivatives: an< 3 The partial 

derivatives of the four constants are determined by the solution of this linear 
system. 


5.4.3 Solution to the Linear System of Unknown Partials 

For the second stage subarc, the solution to the linear system of four 
unknown partial derivatives in the partial derivatives of the four transcendental 
equations is determined by the matrix equation 
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where the equations and ^ arc the same as derived for the analytic 

state partials but are derived with respect to the constant parameters, i.e. 
x = {A/,, C u , C w ). All these terms thus depend on the partial derivatives of 
the common terms a, b, c with respect to the constant parameters. So, 
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Remember that the variation of the terms with respect to the final 

mass is also needed. For the arguments of the inverse hyperbolic sine functions, 

the partial derivatives with respect to the final mass become 
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The partial derivatives of the analytic states with respect to the final mass are 


duj 

dm/ 

dwf 

dm/ 


dh/_ 

dm/ 


C u 1 d^jm/) 

Osfc sj\ + ^2 ( m f) 

9s_ _ Cu/ 1 dS 2 (ro ; ) 

oT +^(m/) dm f 

Ah 1 

o^Tsfc + ^( m/ ) dm s 


{oT? 

Ah 


[m f 


aiaTYfc 

~c w 

a((rT)^/a 
Ah 


, > w ° 
m 0 ) + of 

’ . . / 2 cm/ + b' 

sinh 7 = — 

V v'A , 

/ 2a 4- bm / A 


— sinh 


-l 


sinh 


-l 


— sinh 


-l 


1 


-771/ 


771/\/A 

d3 i(m/) 


{oT) 2 y/c ^1 + Q2( m/ ) dm,/ 


( 2cmo + 6 A 

\T 7 £~) 

f 2 a + bniQ A 
V m 0 \/A / 



51 


_ m 1 dS 2 (m f ) 

1 aT V~a yj\ + %l(m f ) dm f 
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All these relationships are used to determine the coefficient terms of 
the algebraic set of equations. The variations in the constant parameters of 
the zeroth-order two-point boundary value problem with respect to variations 
in the initial states can subsequently be determined. These variations are 
embedded in the quadratures used to calculate the first-order correction to 
the Lagrange multipliers and determine how a change in the initial conditions 
changes the path while flying along a path which will satisfy the terminal 
boundary conditions. 

For the situation where the vehicle has not yet staged, the partial 
derivatives are similar to those shown above but the equations of section 4.4 
which link the two subarcs of the trajectory are used. 



Chapter 6 


Aerodynamic Effect along the Zeroth-Order Trajectory 

Previously the problem of minimizing the fuel into orbit for the flight 
of a rocket in a vacuum over a flat nonrotating Earth was the zeroth-order 
problem, i.e., £ — 0. It was found that this zeroth-order trajectory deviated 
significantly from the optimal trajectory and the resulting correction terms were 
not small as was assumed in deriving the expansion method. To compensate for 
this problem the zeroth-order trajectory needs to be reshaped in order to keep 
the assumed perturbing effects small. One method that might work is to include 
a constraint on the control which will limit the zeroth-order angle-of-attack and 
thus the aerodynamics generated along the zeroth-order path. 1 he problem in 
implementing such a constraint is that the zeroth-order solution must still be 
analytic. Since the analytic solution was found in the local horizon coordinate 
system the control was the pitch angle. Prom the standpoint of the physics 
of the problem, there is no logical constraint which can be imposed on the 
pitch angle. Limiting the angle-of-attack would create a mixed constraint in 
the local horizon coordinate frame involving the state and the control and this 
type of constraint is difficult to solve. A practical and necessary constraint for 
launching a rocket is a dynamic pressure limit, flow such a constraint may be 
incorporated theoretically in the HJB-PDE expansion technique is presented 
in appendix [Cj. But a dynamic pressure constraint arc also does not allow 
an analytic solution to the zeroth-order problem. Therefore, the zeroth-order 
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trajectory was modulated by including aerodynamic terms in the zeroth-order 
problem formulation. This process involved averaging the aerodynamics along 
the vacuum trajectory and solving anew the zeroth-order two-point boundary 
value problem. This technique was suggested by the successive approximation 
method used in [15]. By modelling the aerodynamics as constant terms, closed 
form solutions are still available. This chapter presents the details of includ- 
ing aerodynamic pulse functions averaged in the local horizon and body axes 
coordinate systems. 


6.1 Inclusion of an Aerodynamic Effect in the Zeroth- 
Order Problem 

Instead of assuming flight in a vacuum, the zeroth-order problem is 
now formulated to include aerodynamic terms. Then if e = 0 the equations of 
motion for the zeroth-order problem, valid over both subarcs, become 
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are the assumed lift and drag forces along the zeroth-order trajectory. The 
constant terms A°, A® arc the averaged aerodynamic forces in the x- and 
z-directions. For a vacuum zeroth-order trajectory these terms would be iden- 
tically zero. Nonzero values will be used in order to improve the zeroth-order 
trajectory and keep the perturbation effect due to the neglected aerodynamics 
relatively small compared to the effects due to thrust and gravity. Since these 
terms are added to the zeroth-order dynamics, identical terms of opposite sign 
are included in the perturbation dynamics. Thus their effect is identically zero 
in the full-order system of equations. 


The variational Hamiltonian is altered by the inclusion of these terms, 


e.g. 


T V 

H = — X^V sin 7 4- Av( — cos a — g s sin 7 H ) 

m m 
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V m rn 


(6.3) 


Notice since the pulse functions used in the aerodynamic terms are constants, 
the zeroth-order control law determined by the optimality condition is not 
changed from the solution obtained for vacuum flight. 

(eu) 


Once again the analytic solution to the zeroth-order problem will be 
found in the Cartesian coordinate system. 
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6.1.1 Zeroth-Order Aerodynamic Effect in the Rectangular Coor- 
dinate System 

The equations of motion in a Cartesian coordinate frame become 


h = —w 

T A 0 

w = sin 0 P + g 3 -I - 

m m 

T A 0 

u = — cos 0 p -\ — - (6.5) 

m m 

where the control variable for this problem becomes the pitch attitude 0 V = 
a + 7. The terms A° and A° z represent the constant assumed aerodynamic 
forces along the zeroth-order trajectory in the x- and z-directions, respectively. 
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Figures (6. 2-6. 3) show the aerodynamics averaged over a different number of 
intervals or subarcs. 


The zeroth-order Hamiltonian is 


T a 0 T A 0 
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where A&, A™, and A u are Lagrange multipliers. These Lagrange multipliers are 
propagated by the Euler- Lagrange differential equation A y = —Hj. Thus 


A/t — 0, Au — 0) \pj — A/i 


( 6 . 8 ) 


with boundary conditions 

^h(r/) = Vh, A u (r f ) = i/ U) X w (t/) - u w 
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Fit of Aerodynamic Forces in X-direction 



time 

Figure 6.2: Model for aerodynamic pulses in x-direction 


Fit of Aerodynamic Forces in Z-direction 



Figure 6.3: Model for aerodynamic pulses in z-direction 
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where i/ u , and v w are unknown Lagrange multipliers associated with the 
terminal constraints. Since the aerodynamic effect is added as a constant term 
there is no change in the solution to the Lagrange multipliers or to the control 
from the solution found for a vacuum zeroth-order trajectory. Therefore, the 
zeroth-order analytic state equations become 
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A< = 4a iCi - i> 2 = 4 t-1,2 

and the subscript i refers to the current subarc. More pulse functions could be 
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used to model the aerodynamics in an attempt to capture the effect of the aero- 
dynamics in the closed form solution and thus the path would be broken up into 
smaller subarcs. Note that because the assumed aerodynamics are only con- 
stant terms their effect is an accumulative one. The zeroth-order trajectory is 
altered since the boundary conditions can not be satisfied flying the same path 
as the path flown in a vacuum. The vehicle does not modified its orientation 
instantaneously in order to reduce the aerodynamics that it will encounter, i.e. 
the vehicle cannot predict the aerodynamic effect on the vehicle by its choice 
of angle-of-attack. Thus any change is in the total energy of the system and 
the vehicle is not penalized for flying at large angles-of-attack and for incur- 
ring large drag forces. This can be seen in the new open loop zeroth-order 
trajectory in that the vehicle initially pitches over more than in the vacuum 
solution. But over the entire course of the trajectory the vehicle remains at 
lower angles-of-attack and does not lift up as much in the second stage. If more 
pulses are added the aerodynamics become larger over certain intervals and the 
vehicle reacts accordingly to these regions of large aerodynamic forces. 


6.1.2 FIRST-ORDER CORRECTION TERMS 

The correction terms to the zeroth-order problem can be calculated 
by the quadratures represented in (2.39). Therefore, for the launch problem 
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The first-order term of the optimal return function evaluated along the zeroth- 
order trajectory with initial conditions before staging is written as in (2.26), 
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but separated into two integrals. Only the velocity and flight path angle state 
equations contain the control. Thus, the first-order terms in the expansion of 
the Lagrange multipliers associated with the velocity and flight path angle are 
the only co-state expansion terms needed to construct the first-order correction 
to the zeroth-order control. The partials of P x with respect to the arbitrary 
current conditions, x = (Vo, 70), become 
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Because aerodynamic pulses were added to the zeroth-order dynamics 
the opposite terms are added to the perturbation dynamics such that the over- 
all system equations are unaltered. If the zeroth-order trajectory is the vacuum 
trajectory then the assumed aerodynamic terms ( T> , C) are zero. For nonzero 
assumed aerodynamic forces the new perturbing aerodynamic effect is the dif- 
ference between the actual drag and the assumed drag along the zeroth-order 
path. It is necessary to keep this new perturbing aerodynamic effect small in 
order to accurately approximate the optimal solution. That is the entire reason 
for the inclusion of the aerodynamic pulse functions. The next sections present 
the results for various assumed aerodynamic pulses. 


6.2 Results for the Rectangular Pulse Functions 

It was found that the more pulses used the closer the first-order cor- 
rected solution came to the first-order solution obtained using a vacuum zeroth- 
order trajectory. The best solution for the approximated control was obtained 
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by using one pulse per stage. This seemed to keep the perturbing aerody- 
namic effect small over a larger span of the trajectory. The convergence of 
the Lagrange multipliers up to a first-order approximation using the one pulse 
aerodynamic functions for the zeroth-order problem is demonstrated by the 
plots presented in the Results chapter. Iteration of the zeroth-order trajectory 
for the assumed pulse functions was attempted but it was found that the first- 
order correction terms alternated back and forth between the optimal values 
and the solution based upon the vacuum zeroth-order path. This was a conse- 
quence of the assumed aerodynamics switching between large and small values 
on successive iterations. If large forces were assumed on a particular itera- 
tion than the actual aerodynamic forces along the new zeroth-order trajectory 
would become small and thus on the next iteration the assumed aerodynamic 
pulses would revert to smaller values and therefore the first-order corrections re- 
sembled the solutions obtained using a vacuum zeroth-order path. Attempts to 
average the iterations also proved unsatisfactory. For multiple pulses per stage, 
the averaged iterations did not adequate bring the assumed aerodynamic pulse 
functions closer to the actual forces along the new zeroth-order path. For a 
one pulse per stage solution the iterations could not improve on the solution 
obtained from the first iteration and thus were not worth the computational 
time and effort. In general, assuming more than one pulse per stage and more 
than one iteration caused the first-order corrections to go towards the values 
obtained assuming no aerodynamic forces along the zeroth-order trajectory. In 
a final attempt to lift the vehicle up and keep the vehicle from trying to pitch 
over, aerodynamic pulse functions were modelled as constants in the body-axes 
frame. The next section briefly describes that effort and the results. 
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6.3 Aero Pulses in the Body-Axes Frame 

Because the use of aerodynamic pulses modelled as constant terms in 
the local horizon coordinate system the vehicle did not respond in an instan- 
taneous fashion to the aerodynamics it encountered along a particular flight 
path. To remedy this situation the aerodynamic pulses were modelled as con- 
stant terms in the body-axes frame. Thus there are aerodynamic components 
tangent to and normal to the thrust. Rotation of these forces into the local 
horizon coordinate frame still allows an analytic solution to the zeroth-order 
problem but now the control law becomes a function of the aerodynamic ef- 
fect assumed during a particular interval. This was not the case in using the 
aerodynamic pulses in the local horizon system as presented in the previous 
section. Because of the reliance of the zeroth-order control upon the aerody- 
namic pulses used, the control becomes discontinuous along the zeroth-order 
trajectory. Since the aerodynamic intervals are chosen as functions of a fixed 
time interval the Hamiltonian is also discontinuous across these intervals. The 
integrand used to derive the first-order correction to the Optimal Return Func- 
tion and to the Lagrange multipliers is thus discontinuous and the integration 
of these terms along the zeroth-order path must be broken up according to the 
aerodynamic intervals. The equations of motion in rectangular coordinates for 
body-axes aerodynamic pulses are 
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The terms A ° and A% represent the constant assumed aerodynamic forces 

along the zeroth-order trajectory in the axial and normal body-axes directions, 

respectively. The zeroth-order variational Hamiltonian is 

T 4 - A 0 A° T 4 4° A 0 

H = — XhW+X w ( — sin 0 p +g s cos 0 P )+A U ( -cos 0 P — sintL) 

m m m m 

(6.15) 

The solution for the Lagrange multipliers does not change from the solution 
to the vacuum zeroth-order problem and the multipliers are continuous across 
subarc times, as are the states, since these times are considered fixed. The 
first-order optimality condition produces the following result. 


tan 0 P 


A w(T 4 A 0 ^) + A xj.A°m 
\ W A% - A U (T 4 A^) 


(6.16) 


Using this new control relationship in the state equations the closed form so- 
lution can still be obtained and the states are written as 
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(6.18) 


and the effective thrust f = \J(T + A)\) 2 + (A^) 2 is the magnitude of the sum 
of the thrust and assumed aerodynamic forces. A typical open loop zeroth- 
order trajectory is shown in figure (6.4). While the initial pitch over action 
was curtailed compared to the previous results, the trajectory still deviated 
from the optimal trajectory sharply especially in the regions of high dynamic 
pressure. 


Corrections to the Lagrange Multipliers are made by the familiar 


equation 
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The assumed drag and lift terms are the transformation of the body-axes aero- 
dynamic forces into the wind axes coordinate system, that is, 


Vb = (A 1 ^ cos a — A° n sin a) 

Cb = (A^ sin a + A° N cos a) 


( 6 . 21 ) 
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The correction terms to the Lagrange multipliers based upon the 
zeroth-order trajectory using body-axes aerodynamic pulses did not give any 
improvement over the use of local horizon aerodynamic pulses. If anything 
the solutions obtained were worse since the trajectory was strongly influenced 
(as were the pulse functions) by the regions of high dynamic pressure and thus 
the perturbation aerodynamic effect remained large. The results from iterating 
with the averaged aerodynamic pulses and from averaging the iterations of the 
averaged pulses exhibited the same pattern as the local horizon case. Thus 
one pulse averaged over the first stage came closest to producing agreement 
with the optimal solution. The one positive effect of the body-axes approach 
when used in feedback to generate a trajectory was the elimination of the dis- 
continuities in the control previously found when minimizing the Hamiltonian 
using the first stage aerodynamic model. Unfortunately, the path generated 
did not match as closely the optimal path as the results using the second stage 
aerodynamic model matched. 








Chapter 7 


Results 


In this chapter the approximate optimal solution is compared to an 
optimal solution for the launch of a vehicle in the equatorial plane. While 
previous results for flight in the exoatmospheric regions [16] showed excellent 
matching of the approximate solution with the optimal, problems arose during 
the first stage. First, even at high altitudes where the aerodynamics are indeed 
perturbing effects to the vacuum trajectory, it was found that the linear control 
law derived for the first-order correction to the control (5.1) was in greater error 
than the error in the first-order corrected Lagrange multipliers. As a remedy 
the control was calculated by minimizing the Hamiltonian of the entire system 
using the Lagrange multipliers approximated to first-order. This produced the 
desired effect and the control profile converged to the solution obtained by the 
shooting method. 

The next difficulty encountered was due to the first stage aerody- 
namic model. This model seemed to produce an irregular Hamiltonian. The 
Hamiltonian was badly behaved and exhibited discontinuities in the control 
at various points along the trajectory. The asymmetric configuration for the 
rocket and the cubic spline functions used to fit the aerodynamic data caused 
the Hamiltonian to take on almost identical values for different values of the 
angle-of-attack. This can be seen in figure (7.1) which are plots of the Hamilto- 
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nian versus the angle-of-attack at two consecutive points in the trajectory. The 
sequence shows the Hamiltonian exchanging the location of the minimum be- 
tween positive and negative angles-of-attack. Part of the problem can be seen if 
the drag model is shown for larger angles-of-attack than was presented in chap- 
ter 3.3. Figure (7.2) shows the drag coefficient for different angles-of-attack and 
Mach numbers than would be encountered along the optimal trajectory. Re- 
member the first-order correction terms are based on the aerodynamics along 
the vacuum path but the aerodynamics are not modelled adequately for these 
regions. The drag model of figure (7.2) shows the peculiar nature of the aero- 
dynamics that would be used at the larger angles-of-attack of the zeroth-order 
trajectory. The smooth curve used to model the second stage aerodynamics was 
substituted into the algorithm to eliminate this strange behavior and remove 
the discontinuities in the control. This would prove successful. Figure (7.3) 
compares the drag and lift forces along the first stage of the open loop vacuum 
trajectory using the first and the second stage aerodynamic models. Another 
advantage of using the second stage aerodynamic model can be seen in that 
the drag has been reduced while the lift along the trajectory remains roughly 
the same. 

Overcoming these difficulties still left a problem. The first-order cor- 
rection exhibited a boundary layer type effect near the initial conditions. This 
would occur even if the problem was started at various points in the first stage. 
When the approximation method was used in feedback, this effect would di- 
minish during the trajectory and the solution would converge to the optimal 
solution. In order to eliminate the initial over-corrections of the first-order 
approximation, the zeroth-order problem was reformulated to include an aero- 


Angle-of-Attack a (deg) 


Figure 7.1: Hamiltonian versus Angle-of-Attack at continuous points of the 
first stage 
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Lift Force along Vacuum Trajectory 



Figure 7.3: Comparison of the first stage and second stage aero models along 
the vacuum path 
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Method 

final time 
(sec.) 

final weight 
(lbs.) 

B.C. error 

7 deg 

h ft 

zeroth 

order 

371.50 

322861. 

-0.24 

35. 

first 

order 

369.91 

329293. 

.03 

-.002 

first 

pulse 

369.59 

330576. 

.0001 

.0007 

shooting 

369.57 

330678. 

- 

- 


Table 7.1: Comparison of Results 


dynamic effect. This technique was presented in chapter 6. In this chapter the 
results will be presented along with the results of the zeroth-order solution, the 
first-order solution without the aerodynamic efTect in the zeroth-order problem, 
and the shooting method [17, 18]. 

The trajectories generated by the zeroth-order, the first-order with 
and without zeroth-order aerodynamic pulse functions, and the shooting method 
are shown in figures (7.4-7. 9). Also plotted are the Lagrange multipliers for the 
closed loop trajectory, figures (7.10-7.11). Each technique ran on a IBM 3090 
mainframe computer. Integration was done by an eighth-order Runge-Kutta 
method for the shooting method. The approximate optimal guidance schemes 
employed a fourth-order Runge-Kutta integrator. The approximate method 
used a fixed number of integration steps in the first and second stages with the 
control held fixed over each step. Four hundred steps were used in both the 
first and second stages. The t'ime-to-stage was fixed at 153.54 seconds. 
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All the methods were started at the same initial conditions: to = 35 
sec., ho = 660. ft., Vo = 9406. ft/s, 70 = 58. deg., mo = 3021 107. AAibs., 
0 o = -79.0 deg., and x = <t> = 0.0 degrees. The terminal constraints to 
be satisfied are hj = 486080. ft., V } = 25770. ft/s, and 7/ = 0.0 degrees. 
The results are compared in Table (7.1). The solution shows the approximate 
optimal guidance law using the first-order correction term matches the control 
and state trajectories of the shooting method. Initially only the first-order 
correction with the aerodynamic pulse generates a nearly optimal trajectory. 
The cost obtained by the two techniques is nearly identical. The final weight 
using the shooting method was 330678. lbs. at a final time of 369.57 seconds. 
The final weight was 330576. lbs at a final time of 369.59 seconds when using the 
first-order approximation. The zeroth-order solution shows a greater variation 
in the control from the optimal control. The final weight obtained was 322861. 
lbs. at a final time of 371.5 seconds. The zeroth-order solution also does 
not satisfy all the boundary conditions as closely as the optimal and first- 
order solutions, with an error in the final flight path angle of -.24 degrees 
and an error in the final altitude of -1-40 feet. Because of this error in the 
terminal constraints, large angles-of-attack can be seen in fig. (7.4) for the 
zeroth-order solution in attempting to meet the terminal constraints. The first- 
order correction picked up most of the deviation of the zeroth-order trajectory 
from the optimal trajectory and as a result the boundary conditions are met 
more closely with a better behaved control. The most important aspect in 
obtaining good results is the convergence of the Lagrange multipliers to the 
optimal Lagrange multipliers. With the use of the aerodynamic pulses the 
flight path angle Lagrange multiplier approximated to first-order shows good 
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time 


Figure 7.4: Angle-Of-Attack vs. Time 
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Figure 7.6: Altitude vs. Time 
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Figure 7.7: Velocity vs. Time 
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Figure 7.8: Flight Path Angle vs. Time 
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Figure 7.9: Dynamic Pressure vs. Time 






Figure 7.11: Flight Path Lagrange Multiplier vs. Time 
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Method 

zeroth 

order 

first 

vacuum 

first 

pulse 

shooting 

CPU time 
(sec) 

49. 

304. 

344. 

426. 


Table 7.2: Comparison of computation time 


agreement with the optimal solution. A last point about these result is that the 
inclusion of the rotation of the Earth in the problem is expected to continue to 
reduce the time of flight and consequently increase the final weight available at 
orbital insertion. 

The convergence of the asymptotic expansion is indicated by the re- 
sult of the first-order solution in comparison with the shooting method so- 
lution, thereby precluding the need to include higher-order correction terms. 
This convergence is tentative since it took the inclusion of the aerodynamic 
pulse functions in the zeroth-order problem to achieve the best results. Alas 
the convergence properties when using these pulses cannot be guaranteed or 
even quantified. Finally, since this algorithm is being proposed as a real-time 
guidance scheme the computational time that was needed to generate the entire 
trajectory by each method is presented in Table 7.2. While none of the codes 
have been optimized for computational efficiency, the use of quadratures does 
decrease the time needed to solve the launch problem in comparison to the 
shooting method. It should be noted that the flight time is approximately the 
same as the cpu time for the first-order approximation methods and that the 
shooting method was given a good initial guess (nearly converged) of the un- 
knowns. As expected, the zeroth-order analytic solution was found extremely 
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quickly. The introduction of the aerodynamic pulse functions into the method 
caused a modest increase in the computation time. 



Chapter 8 


The Relationship between Calculus of Variations and 

the HJB equation 


In this chapter the relationship between the Hamilton-Jacobi-Bellman 
approach and the calculus of variations approach is presented. First, the HJB 
expansion method is described in more detail in order to explicitly show the link 
between it and the perturbation of the canonical form of the Euler- Lagrange 
equations. The similarity of the terms involved and of the two solution tech- 
niques is shown. Next the solution of linear, first-order, partial differential 
equations is described. The significant result of the solution process is that 
the solution to a partial differential equation is equivalent to the solution of 
the characteristic curves represented by a set of ordinary differential equations. 
A simple derivation of the Lagrange multiplier differential equation from the 
HJB-PDE is also included. Lastly, the formulation of the ALS problem along 
with the results obtained when using the expansion of the calculus of variations 
method are presented. 

8.1 Correction Terms to the Lagrange Multipliers 

In Chapter 2 the equation for the Lagrange Multipliers (the change 
in the cost due to a change in the initial state) was determined to be 
dPi(x,t) f T fdRi(y,r,P 0l ) dr f 

— “ -y, — dx — iT ~ R '\’’S 
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where the terms axe evaluated along the zeroth-order path and the higher-order 
expansion terms Pi(y/, tj) = 0 (z = 1,2, . . .) at the terminal boundary. 

Now recall that the integrand terms were 


P\(y,r,Po x ) = -Pj x (g(y,u,T) + fi(y,u,r)) 

dfli( 2 /,T,P 0l ) _ t) + fi(y,u,T))] 

dx dx 


( 8 . 1 ) 

( 8 . 2 ) 


where the expansion term in the primary dynamics f\ is 

df(y,u,T)\ 


fi(y,u,r) = 


du 


u i 


£=0 


(8.3) 


Along the zeroth-order path the optimality condition Po x f u — 0 eliminates the 
fi term in the integrand R\ 


Ri(y,r,P bj = - Po x g(y,u,T ) 


is 


Therefore, the first-order term in the Lagrange multiplier expansion 


dPi(x,£) 

dx 



dg^dy dg_du 
dy dx du dx 


dr 


(8.4) 


Now the equation for the zeroth-order control can be written as a 
function of the independent variable, of the states, and of the partial derivative 
of the zeroth-order Optimal Return Function with respect to the states. 
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du 

dx 

du 

dx 


Thus the variation of the control with respect to the initial state is 
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Substituting this equation into the first-order Lagrange multiplier 
equation results in 
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(8.5) 


Note the notation used here is that the partial derivatives are taken with respect 
to an individual initial state x not the initial state vector. 


The integrations thus depend on the variation of the zeroth-order 
states and the zeroth-order Lagrange multipliers due to variations in the initial 
state x, i.e. the terms and represent the state transition matrix. This 
matrix can be obtained from the zeroth-order analytic solution. 
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where the equation for the control variation has been used. Similarly for the 
Lagrange Multipliers, 
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These coupled equations could be written in matrix notation as 
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where A yy , A y p 0i , Ap 0xV , Ap 0i p 0z are the coefficients of the differential equa- 
tions presented above. 

The change in the parameters associated with the zeroth-order two- 
point boundary value problem due to a change in the initial states can be 
determined by the variation with respect to the initial states of the terminal 
boundary conditions and of the Hamiltonian at the final time. For every change 
in the initial state the transversality conditions must still be satisfied. So, 
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drj dx 


dP 0l dy + dP 0x dTf 


dy dx dr/ dx J 
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= ^0 ( dfar, T d ^r, \dy ( d< k>'t r gjs \ djy 

°v dx \ dy 0 dy ) dx \ dr/ 0 Otj ) Ox 

The variation in the parameters associated with the two-point bound- 
ary value problem with respect to variations in the initial states can be found 
from these sets of equations. Notice these equations are linear in the unknown 
parameters f^, and 

8.2 Expansion of the Euler-Lagrange Equations 

This section attempts to relate the results of the expansion of the 
HJB-PDE to the results derived from expanding the ordinary differential Euler- 
Lagrange equations as was done in [19]. 

The states, control, and Lagrange multipliers are all expanded in an 
asymptotic series with expansion parameter e. Thus, 


y(r) 

= yo(x) + £ 2 /i(t) + 0{e 2 ) 

(8.6) 

u{t) 

— uq(t) + eu\(r) + 0(e 2 ) 

(8.7) 

A y (r) 

— A 0v (t) + eAi y (r) + 0(e 2 ) 

(8.8) 


These expansion equations are used in a Taylor series expansion of 
the dynamics: 

Jo = f(y,u,r) (8.9) 

/i = fu{y,u,r)u^ + f y (y,u ) T)yi (8.10) 


go = g(y,u,r) 


( 8 . 11 ) 
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9\ = 9u(y,u,T)ui + <] y (y,u,T)yi (8.12) 

where the notation is the same as previously presented, i.e., the state and 
dynamics y, f, and g are n-vectors and the control u and independent variable 
r are scalars. 

When expanded in the small parameter e , the optimality condition 
\ T (df(y,u,T) , dg(y,u,r)\ 

A * ( dn +e ~~ au ) =0 

becomes 

df(y,u,T) 

\T df(y,u,r) xT ["d/i(y,u,T) dy(y,u,r) 

A, » _ " au + M Tu — + -~du J = 0 

after the coefficient terms of like powers of e have been collected. 

Thus the first-order expansion term for the optimal control is 

u i = ~ [A£/uu (y.ti.r)] ‘ [lu(y,u,r)\ ly + Xlj uy (y,u,T) yi + \£ v g u (y, u, r)] 

This result is substituted into the Euler-Lagrange ordinary differential equa- 
tions. 


8.2.1 Expansion of the State Equations 
The differential state equations 

$ = /(?/> u > T ) + £ 9(y . w, r) 
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with initial conditons y{r = t) = x, can be expanded in the small parameter e 
as 

yo + eyi + e 2 ij 2 + ... = fo + s{fi + go) + ^(fi + Q\) + . . • 

Collecting like powers of e leads to 


f(y,u,r) 

(8.13) 

U(y,u,T)u i + f y (y,u,T)yi 


+g(y,u,r) 

(8.14) 


The initial conditions are yo(t) = Xq , V\{t) = x\ = 0 . 

Thus, 

£o(r) = f(y,u,r) 

2/i( r ) — | fy ~ fu (^ 0 ,/uu) \) v /uy 2/i — /u (^o y /uti) / u A Iy 

d" 9 ~ fu ^A 0y /uu) A q v 5u | 

when the expansion terms of the control law are inserted into the state equa- 
tions. 

8.2.2 Expansion of the Lagrange Multiplier Equations 

Expanding the differential equation for the Lagrange Multipliers 

A u = —(fy + e 9^) Ay 
produces the following set of equations 
Ao v = -/y T ( 2 /.u,r) A 0v 


(8.15) 
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A 'v = ~fy(y> U > T )^l v - fT v \ - 9y(y,V-,T)\o v 


(8.16) 


Therefore, 


Ao„ = - 


d/(y,«,T) ' 

dy 


Ao„ 


• _ dfjy 1 u 1 r) T 

ly dy A|y 


d fT (y. u, t) + dy T (y, u, t) 




dy 


(8.17) 
A 0y (8.18) 


Consequently, the differential equations which describe the first-order expan- 
sion terms are 


d 

dr 


[ yt ' 


A;y 

St 


!/. 1 

1 

s> 

t 


^ v y 

^A„A y 


>* 

< 

1 


where 


Ayy — I fy fu (^ v /uu) ^0 y /t 

A/A„ = -/u(AJ v A u ) _ 1 /u T 


l A v y ~ 




yfyy \ y fuy [\l v U) \lj y 

/ s r-AoVu V (A 0 v uu )- i /j 


(8.19) 

( 8 . 20 ) 
( 8 . 21 ) 

( 8 . 22 ) 


The solution to this set of coupled differential equations is 

. w ] = [ l\% ] + C ^r,r)G(r)d, 

with G{t ) representing the forcing terms. 


G(r) = 


— A 0y + A 0j( /uy (Ao y /uu) 9u Ao y 
This is the same transition matrix as derived for the variations of the Optimal 
Return Function. Notice that the forcing terms are also the same as those 
derived for the variation of the Optimal Return Function. 
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8.3 Expansion of the Boundary Conditions 

In order to expand the terminal boundary conditions properly the 
states must be expanded about the zeroth-order path and the zeroth-order 
final time. Thus, 


7/ — tq { + et\ { + 0{e 2 ) (8.23) 

y( T f) = yo(r Qf ) + e [ 2/1 (r 0/ ) + 2 /o(t 0/ )ti / ] + 0{e 2 ) (8.24) 

= 2 / 0 / + £y\ f 4- 0{e 2 ) (8.25) 


Next expand the terminal boundary conditions with respect to the 
small parameter e: 


v(y/> T /) = = 0 

t=0 

= [^0 (y/,T f ) + e'H\(yj,T f ) + £ 2 ^ 2 (y f ,T f ) + 0(e 3 )] 

The Taylor series expansion of the terminal boundary conditions is 

1 d X {ifj f Tj ) N / 00 W / 1 'I-O^ \ / 00 


^(y/. r /) = H 


i= 0 


(r 


dy'r 




= *(yf> T f) 


-\-£ 




=0 dy f 


e=0 


,J, ' +e d7- 


dV 
7 


£=0 


T\ f + 0(e 2 ) 


Collecting like terms in the expansion parameter e results in 


Vo{yf,T f ) = ty(y f , T f ) 


= 0 




dy f 


E —0 

dV 

2/1/ + o — 

£ =0 OTf 


(8.26) 


£=0 


T\ f = 0 


(8.27) 
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8.3.1 Expansion of the Transvcrsality Conditions 

Expanding the Hamiltonian A ^ (/ +eg ) in terms of the small parameter 
e produces the following equations 

Hamo = Aj 'j{y,u,r) 

Ham i = \o v gi.y,u,r) + \^J{y,u,T) + XoJ y (y,u,T)yi (8.28) 


with boundary conditions at the terminal time 


IIamo(ro f ) 

Ham\(rf) 


d(p(yQ r T 0f ) r ( 9 ^(y 0 /I r 0/ ) 

V n 


dr f 


drt 


iyof i tqj ) d(j) T f (2/0, > ) 7’ d't (2/0, > to, ) 


dr, 


-2/i 


dr/ Tl/ 1/1 dr; 


f d^j/, (z/o y ) ^0, ) 7’d'I , r, (2/0, > ^C),) 

— a^ m > ~ 1/0 a^ n < 


Equating the first-order expansion terms of these relationships, the 
equation for the first-order Hamiltonian at the terminal boundary becomes 


Ao y ( r O/)5(2/o/> lt o/) r o / ) + A£(t 0/ ) + Ao, (t 0/ )t 1/ /( 2/0,, u 0/) r 0/ ) 

1 \T / \ e / ^ ^V/(2/0/i r 0 r ) dfir/iyOfiTOf) 

+K y (To f )f y (yo / ,uo f ,T 0f )y lf = — 1 — — 


2/ir - 


&Tf { 


T d^ yf (yo r T Qf ) _ r 9^ T/ (7/ 0/ ,Tb / )_ _ T 5^(2/ 0/ ,r 0/ ) 

— i/ n ^ ?/l r 7 : T \ f “ 


"° dr/ yi/ ^ drf 


Now proceeding in a similar fashion, the terminal Lagrange Multi- 
pliers, determined by the relationship A^(r/) = (p y/ (yf,Tf) + j/ T 'I' !(/ (y; , r;) , 
become 

30(2/0 ,, t 0/ ) _ r d v L(y 0/ , r 0/ ) 
dyf U ° dyj 


Ai„(to,) 


(8.29) 
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, , , , x - , N 9(j) y (y 0 t q ) d<p T (y 0 ,t 0 ) 

Xo v (r 0f ) + X 0y (r 0f )T lf = L—L—^ y Xf dy ; 


■ dVjyo^ro/) 


r d^ yf (yo f ,r Qf ) n 

~ u ° — W, — Vh 


^d^r f (yo f ,T 0f ) 


The expansion of the t ran s vers ality conditions and the terminal bound- 
ary conditions produces the boundary conditions for the first-order two-point 
boundary value problem. As in the case of the variation of the Optimal Return 
Function, the unknown parameters, (yi f , v\, and T\ f ), can be found by solving 
a set of linear algebraic equations. 


8.4 Solution to the First-Order Problem 

The solution to the first-order two-point boundary value problem is 


found by use of 


] = ' 0 [ l% ] + r *^0„T)G(T)dl 


(8.30) 


subject to the terminal boundary conditions given for y\{rQ f ) and Ai B (ro / ) in 
Eqs. (8.27 - 8.30) and subject to the Hamiltonian transversality condition 
Eq. (8.29). Also recall that the initial states are considered known and are 
zeroth-order terms. Thus the first-order initial states are set equal to zero, 
y\{t) = 0. The unknowns which need to be found are the initial and terminal 
first-order Lagrange multiplier terms, i.e. Ai v (f) and V \ , and the first-order 
term in the expansion of the final time r t/ . 

For a trajectory that includes the staging condition, the terminal 
conditions remain the same but the form of the solution becomes 

yi( T o f ) 1 _ * / T n* .x T yi(0 

Ai v (r 0/ ) j ~ [ A ly (0 
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/■^Oy fTatagc 

+ / $a 2 (to t)G(t)cLt 4- ^a 2 (t 0 T stage ) / ^ A x {r sla g e ,T)G(r)dT 

Jr.tagc Jt 

where the state transition matrices and ( I> 4 2 represent the state transition 
matrix over the first stage and second stage subarcs, respectively. 

Notice that the form of the solution determined by expanding the 
Euler- Lagrange equations and the terminal boundary conditions in terms of 
the states, the control, and the Lagrange multipliers is equivalent to the solu- 
tion found by the expansion of the Hamilton- Jacobi-Bellman partial differential 
equation. The forcing terms and the transition matrix used in the quadratures 
are the same. The first-order boundary conditions derived in this section are 
identical to the variation of the zeroth-order boundary conditions which are 
used in the HJB expansion method to determine the change in the parameters 
of the zeroth-order solution with respect to a change in the initial states. In 
the HJB expansion the variations and arc dependent on the changes in 
the constant parameters or the constants of the motion. Any admissible vari- 
ation in the initial conditions must still generate a trajectory which satisfies 
the terminal conditions. Thus the variations in the boundary conditions with 
respect to changes in the initial states determine the change in the constants 
of the motion. And these changes are embedded in the solution of the vari- 
ations in the state and Lagrange multipliers, ^ and which are used to 

generate the first-order correction terms. In contrast, the first-order boundary 
conditions derived from the Euler- Lagrange equations (which are equivalent to 
the variation of the zeroth-order boundary conditions) are explicitly used in 
solving the two-point boundary- value problem (8.30). 

The similarity of the two techniques is not surprising since solving the 
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HJB-PDE is equivalent to solving the first-order Euler- Lagrange equations. 
The reason for this equivalency will be presented more clearly in the next 
section. 

8.5 Solutions to First-Order Linear Partial Differential 
Equations 

The solution to first-order partial differential equations is described 
in [11, 20]. In this section it will be shown that the canonical system of Euler- 
Lagrange differential equations is identical to the system of characteristic ordi- 
nary differential equations used to solve the partial differential equation. This 
is presented for the partial differential equation in two independent variables 
but the case of n-independent variables is a straightforward extension. First 
consider a partial differential equation of the form 

F(x, t, P x , P t ) = a(x, t)P t + b(x, t)P x - c(x, t) = 0 (8.31) 

where a, b, and c are given functions and arc considered continuous, as are 
their first derivatives, in the region of interest. The solution of this partial 
differential equation is called the integral surface and is denoted as P(x,t). 
Since the coefficient terms ( a,b,c ) are not explicitly dependent on the solution 
P this is a linear rather than quasi-linear partial differential equation. The 
tangent plane to the integral surface P at the point Q(x, t, P) is defined by 
the relationship (8.31) and the normal to the tangent plane is given by the 
directions P x , P t , and —1. The partial differential equation implies that the 
normal to the integral surface < P x , P t , —1 > is perpendicular to a vector 
< a, b, c > and so the vector < a, b, c > must be tangent to any integral 
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surface at point Q. The vector < a,b,c > is called the Monge vector. The 
tangent planes to all integral surfaces through the point Q{x,t , P ) belong to a 
family of planes which are described by the relation for the normal as 

dt : dx : dP = a : b : c 


The direction of the Monge vector at the point Q forms a characteristic line 
element ds. The directions of the Monge vectors form a directional field in 
the (x,t, P)- space. To solve the partial differential equation (8.31) the surfaces 
which fit the Monge vector must be found. Every surface whose tangent plane 
is tangent to the Monge vector at the point is a solution to the partial differen- 
tial equation. The characteristic curves of the partial differential equation are 
the integral curves of the direction field and are defined by a set of ordinary 
differential equations. If the characteristic curves are considered a function of 
a parameter s then along the curves the characteristic equations become 


dt 


— a 


dx . 

Ts =b 


dP 


— C 


(8.32) 


ds ” ds ds 

Thus a general solution surface can be generated independently of the initial 
data as a one-parameter family of characteristic curves. 


For the initial value problem the manifold of possible integral surfaces 
can be created and the unique solution depends on the initial conditions of the 
problem, i.e. y(s = 0) = x, r(s = 0) = t, P(x, t) =constant. Starting with a 
curve S in space the solution to the partial differential equation is sought. The 
curve S is projected onto the (x, f)-space and an integral solution P(x,t) is to 
be found, see figure (8.5). Through each point of the space curve a family of 
characteristic curves can be generated according to the characteristic equations. 
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These curves form a surface and all characteristic curves lie within this integral 
surface. 

The solution is sought at a point off of the initial data curve and in 
the direction of the characteristic curve. Thus the solution becomes a function 
of two variables, the initial state (y(s = 0) = x, t(s = 0) = t) and the running 
parameter s along the characteristic curve. As a consequence the solution to the 
integral surface P(x, t ) can be written as a function of s only with x, t replaced 
by their respective solutions along the characteristic direction and with fixed 
initial conditions at s = 0. This can be done if the characteristic solutions for 
x, t can be inverted to obtain functions dependent on y(s = 0) — x, t(s = 0) = 
t. The transformation between the two sets can only occur if the Jacobian is 
nonzero. In this case a unique solution exists to the initial value problem. 

If the characteristic curve and the projection onto the x,t plane of 
the tangent to the curve S are identical then the curve S is a characteristic 
curve. Mathematically this happens if the Jacobian is equal to zero for every 
point along the curve. This is the relationship that was obtained in chapter 
2 for the Hamilton- Jacobi- Bellman equation where a = 1 and b = f(x,u,t). 
The solution obtained along the characteristic by definition must stay along the 
initial data curve. The implications of this result are that an infinite number 
of solutions exist for the integral surfaces which solve the partial differential 
equation and which pass through the curve S. Since the ordinary differential 
equations for the characteristic curves require the integration of a,b, and c 
which are known data along the the space curve S and since the projection of 
S in the x, t plane is the curve S itself, the integral solution P(x, t) can be found 
by integration °f — c. The unique solution to the problem is determined 
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by the initial conditions of the ordinary differential equations and the terminal 
boundary conditions such that P can not be specified arbitrarily at the initial 
point y(s = 0) = x, r(s = 0) = t. 

In [10] the relationship to the calculus of variations approach is simply 
derived. If a particular trajectory is the optimal trajectory then the Hamilton- 
Jacobi-Bellman equation must hold at each point and thus can be written with 
respect to the states y and independent variable r. The fundamental equation 
is 


Py (y, r) (f{y, u, t) + eg(y, u, r)) + P T (y, r) = 0 

and the optimal return function depends only on the cost at the terminal 
manifold. Therefore, the optimal return function is constant and the total 
time derivative must be zero at each point in the path. The partial derivative 
with respect to y is 


o = (l T (y,u,r)+eg T ( v ,u,T))^ + ^- + P^ 

■ pT ( jjg 9g(y,u,r) du\ 
y \ du dy du dy ) 

By the chain rule for differentiation 


(df(y,u,T) dg(y,u,r)\ 

V dy +£ dy ) 


d ^y ( tT/ x , T f ,\dP v dP T 

Thus a system of ordinary differential equations is obtained which is evaluated 
along an optimal trajectory and is satisfied by the partial derivatives of the 
Optimal Return Function P(y,r). 


iP v _ _ pT ( df {y , u,t) , dg(y,u,T)\ 

dr f dy +£ ay ) 

p t ( df(y,u,r) du dg(y,u,T) 

y V d u du dy) 


(8.33) 
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Now the optimality condition is 

nT f d f(y> u > T ) , .9g(y,u,r)\ _ n 
^ \ du du ) 

Therefore, 

dp y _ _ pT ( ^/(y- u » r ) , dg(y,tt.r) \ (8 34) 

rfT" '» V d 'J d y ) 

which is the familiar Lagrange multiplier rule for the optimal trajectory. 

The material presented is intended to clarify the relationship between 
the solution to the first-order, linear, partial differential equation which is 
the Hamilton- Jacobi-Bellman equation and the solution to the Isuler-Lagrangc 
equations. The point to remember is that the solution to the partial differential 
equation is given by characteristics generated by ordinary differential equations 
which are the equivalent to the canonical Euler- Lagrange equations. 


8.6 Formulation of First-Order Correction Terms for 
the ALS Problem 


In this section the solution for the ALS problem using the new per- 
turbation method , i.e. the perturbed Euler-Lagrange equations, is presented. 

The state transition matrix is determined by integrating 


ar 


^ yy -^yAy 

h„y -^AyAy 


( h(T,0. 


where the A matrix was presented in Eq. (8.22) 


<L(M) = / 


For the ALS problem the primary and perturbation dynamics are 



’ w ' 


■ g s _ 2W sin 0 p ■ 

II 

?r 

u 

— 

^ cos Or) 

m r 


_ h . 


— w 
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and 

( 2 »+; 
r e +/i) 
_ _n 
rn 


where and A z represent the aerodynamic forces in the x- and z-directions, 
respectively. Also n is the number of engines per stage and p and A e are the 
pressure and the engine exit area. Remember that the thrust was modelled as 
T = T vac - npA e . 




The first-order partial derivatives of the primary dynamics are 

0 0 0 — cos 0 P 

f y = 0 0 0 / u = sin 0 P 

_ -1 0 0 J [ 0 

and all the second-order derivatives are identically zero except / uu , 

r ~T sin 

^ cos 0 P 
0 

Therefore the matrix A in the differential equation defining the state transition 
matrix becomes 

0 0 0 Q14 U]5 0 

0 0 0 (Z24 &25 0 

.-1000 00 

0 0 0 0 0 1 

0 0 0 0 0 0 

.0000 00 . 

with 



Tvac cos 2 Op 

m (A u cos Op — \ w sin 0 P ) 

°15 = &24 = 


_ Tvac sin 2 Op 

m (A u cos 0 P — X w sin 0 P ) 
T vac cos Op sin 0 P 
m (A u cos Op — A^ sin 0 P ) 
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An analytic form for the state transition matrix can thus be obtained. 

0 0 <f > 14 

1 0 24 

<*>34 

1 


4 >(r ,0 = 


and 


<f>15 = $24 — 


1 

0 

-(r — t) 0 1 

0 0 0 

0 0 0 

0 0 0 


0 

0 


$15 

$25 

$35 

0 

1 

0 


$16 

$26 

$36 

(r - t) 
0 
1 


$u = -- 


ci _ 

go ?! 2 

_2 Cl 

oa A 

C U C W 

ama ?/ 2 

2 C U C, 


sinh 


-l 


2 a 4 - bm ' 
V mV A , 
( bem - 2 ac 4 - b 2 ) 


— sinh 


-l 


2ci 4" brriQ 


y m 0 V A 
( hem, — 2 ac 4 - 6 2 ) 
\/ cm 2 + bm 4 - a y/cm 2 4 - 6 m, 4 - a 


sinh 


-l 


2 a 4 - brn 


+ 


$16 — 


ama A 

2 C U A„ 
ark A 

Cfoi 

ornd s l 2 

2 C 2 m ^ 
ama A 

2 C 2 


mV A 
(tem — 2 ac 4 - 6 2 ) 
Vcm^TlOTrfa 

( 2 cm 4 - 6 ) 

\/ cm 2 4 - Inn 4 - a 

sinh -1 


sinh 1 
( 6 cm, 


" 2 a 4 - 6 mc 
v rn oVA 
- 2 ac 4 - 6 2 ) 


2 a 4 - 6 m \ 


\Jcrn 2 4 - 6 m, 4 - a J 

(2 cm, 4 - 6 ) 

\Jcrn 2 4 - 6 m, 4 - a_ 

2 a 4 - 6 m 0 ' 


m 


v/a ; 


sinh" 


$25 — ~ 


+ 


am A 

c 2 

aa 3 / 2 

2 Cl 

aa A 
26 


aA 


\/ cm 2 4 - 6 m 4 - a 

( 2 cm 4 - 6 ) 

\/ cm 2 4 - 6 m 4 - a 

, _i / 2 a 4 - 6 m \ 
sinh 7 =- 

V mVA J 

( ban — 2 ac 4 - 6 2 ) 

\/ cm 2 4 - bm 4 - a 

( 2 cm 4 - 6 ) 


moVA 

(bem — 2 ac 4 - 6 2 ) ( 6 cmi — 2 ac 4 - 6 2 ) 

cm 2 4 - 6 m, 4 - a 

( 2 cm* 4 - 6 ) 

\jcm 2 4 - 6 m, 4 - a_ 

/ 2 a 4 - 6 m 0 \ 
V m 0 \/ A y 

2 ac 4 - 6 2 ) 


- sinh 1 
( 6 cm, — 


y cm 2 + brrit 4 - a J 
( 2 cm, 4 - 6 ) 


V cm 2 4 - 6 m 4 - a y / cmf+ 6 m 7 +aj 
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2c 


(2a 4- brri) 


(2a 4- 6m,) 


aA [Vera 2 + bm + a ^cm 2 + 6m, + aj 


a c m m, 
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2C u C w m x 


sinh 


-l 
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VA 


- sinh 
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'2a 4- brriQ 
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amaA 
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+ 
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20A* 

crm 2 A 

c 2 
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2Clb 


s~-t . 

Cu,. 4- —rn x 
rn 
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amaA [ 
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amaA 
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vX? 4- bm 4- a ^cm 2 4- 6m, + a J 
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m sinh =- - sinh =- 
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ama 


3/2 
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sinh 
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The first-order partial derivatives of the perturbation dynamics are 


9y = 


' 

Tf. dA z 

Is. 

< 1 DA Z _ 

2u \ 

dw 


mh 3 dw 

h. 

t m tfu 

(rc+k)J 

dh. 

Ls. 

( 1 dA x | \ 

Is. 

1 9A X , 

w \ 

du 

h M 1 

\ m dw ’ {r c +h)J 

h, 1 

\m du ' 

(r e -f/i)y 

dh 

_ 

0 


0 


0 


and 


9o„ = 


mh , 


'l£%, + n vAe cos Q p ; 
-gf- + npA e sin 0 P j 
0 


where 


dw 

dh 

du 

dh 


1 dA, 


4- 


u 


2g,rl npA' s . n 


m dh (r e 4- h) 2 (r e -i- h) 2 mh 


1 dA 3 


f e I 1 us\x UW ups- i e 

Y 1 — — - /- , + ~HT~ cos °v 


npA e 


(8.35) 

(8.36) 


jn .dh (r e 4- h) 2 mh v 

The forcing terms in the quadratures pertaining to the first-order 
correction terms can now be calculated. Recall 

5-/«( a o v /uu) 1 glK y 

~9y ^ Oy + ^o y /«y (^o y /uu) A 0 , 


G(r) = 


which becomes 


G(r) = 


9-fu(^oJu U ) 1 9u^0 


~9yK 


The aerodynamic forces are considered positive in the x- and z-directions 
respectively, and are 


A x = — 0.5p5(u 2 4- w 2 ) [C d (M, a:) cos7 -I- (^(M, a) sin 7] (8.37) 

A z = 0.5pS(u 2 4- w 2 ) [Cq(M, a) sin 7 — Ci(M, a) cos 7] (8.38) 

The angle-of-attack is a function of the pitch angle and the flight path angle, 
01 — 9 P — 7. The atmospheric density and pressure are modelled as exponen- 
tials, i.e. p = p s exp(—h/h s ) and p = p a exp(—h/h p ). The mach number is a 


105 


function of velocity and altitude with M — ( u 2 + w 2 )/sos. The specd-of-sound 
is calculated by sos = ^F^, where F is the specific heat ratio for air and is 
assigned a value of 1 .4 . Lastly, the flight path angle is represented in the Carte- 
sian coordinate system as tan 7 = — From these relationships all the partial 
derivatives needed to calculate the forcing terms G(r) can be determined. 

8.7 Results 

The solution to the launch problem was first attempted for initial 
conditions associated wit!) staging. At these altitudes the aerodynamic forces 
are small enough that they may correctly be consider perturbation terms. The 
results of the new peturbation method (expansion of the Ruler- Lagrange equa- 
tions) show excellent agreement with the optimal solution. Note that the entire 
first-order correction is available since this method is valid as an open loop solu- 
tion as is shown in figure (8.2). The results also agree exactly at the initial point 
of the path with the previous results using the old method (HJB). Table (8.1) 
lists the rclevent values. 

Next, the solution for initial conditions at a time of 35 seconds was 
sought. Once again the solution via the new method matched exactly the re- 
sult obtained using the old method. To obtain agreement with the optimal 
trajectory, the aerodynamic pulse functions were utilized in the same manner 
as previously discussed. Consequently, the first-order solution closely approxi- 
mated the optimal solution. The values at the initial point are also included in 
Table (8.1) and plots of the open loop profiles are in figure (8.3). The solution 
in a feedback configuration is presented in figures (8.4-8. 6). Presented are the 



Method 

t 0 = 35. 

To = 153.54 

Pv 

P 7 

Pv 

Py 

new 

first 

-2.2535 

-843.12 

-0.8547 

-908.35 

new 

pulse 

-1.3925 

-153.50 

“ 


HJB 

first 

-2.2535 

-843.12 

-0.8547 

-908.36 

HJB 

pulse 

-1.3954 

-153.82 


“ 

shooting 

-1.2752 

-139.48 

-0.8151 

-860.63 


Table 8.1: Comparison of open loop results 


Method 

final time 
(sec.) 

final weight 
(lbs.) 

B.C. error 

7 deg 

h ft 

new 

first 

369.91 

329295. 

0.0026 

.219 

new 

pulse 

369.59 

330578. 

.0014 

-.144 

HJB 

first 

369.91 

329293. 

.03 

-.002 

HJB 

pulse 

369.59 

330576. 

.0001 

.0007 

shooting 

369.57 

330678. 

- 

- 


Table 8.2: Comparison of closed loop results 
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profiles of the flight path angle Lagrange multipliers, the velocity Lagrange 
multipliers, and the control over the entire trajectory. In the closed loop solu- 
tion the values are practically identical with a slight difference in the accuracy 
upon which the terminal conditions are met. In order to be consistent with the 
results chapter, the second stage aerodynamic model was used throughout the 
flight. Table (8.2) verifies that the results are the same for the two perturbation 
methods. 

To summarize, the expansion of the Hamilton-Jacobi-Bellman equa- 
tion is equivalent to the expansion of the Euler-Lagrange canonical equations 
with respect to the states, control, and the Lagrange multipliers. The theoreti- 
cal and geometric concepts behind the solution of first-order partial differential 
equations were also discussed. The reason for the equivalency of the two meth- 
ods is that the result obtained by solving the partial differential equation and 
differentiating with respect to the initial states is identical to the result ob- 
tained by the solution of the ordinary differential equations that represent the 
characteristic equations for the partial differential equation. For the Hamilton- 
Jacobi-Bellman equation the characteristic equations are the Euler-Lagrange 
equations. As expected, the solutions obtained using the two perturbation tech- 
niques are identical. While the calculus of variations approach took a longer 
amount of computation time, the entire open loop trajectory can be generated. 
Because of this fact the update to the feedback solution need not be com- 
puted as often and thus the overall computational time can be reduced. At the 
expense of this speed comes the additional burden of integrating a state tran- 
sition matrix rather than calculating the partial derivatives of the zeroth-order 
solution. While the state transition matrix approach is easier to understand 
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time 


Figure 8.3: Open loop solution for Lagrange multipliers at first stage initial 
conditions 
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Figure 8.4: Closed loop solution for flight path angle Lagrange multipliers 
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Velocity Lagrange Multiplier 



time 


Figure 8.5: Closed loop solution for velocity Lagrange multipliers 
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than the embedded nature of the HJB solution, the state transition matrix can 
be more difficult to obtain than the corresponding partial derivatives needed 
by the HJB expansion method. Also note that the Hamilton-Jacobi-Bellman 
equation can be written as a stochastic equation and can thus handle random 
disturbances. 



Chapter 9 


Conclusions 


The technique for applying the expansion of the Hamilton- Jacob i- 
Bellman partial differential equation to derive a real-time guidance scheme has 
been presented. The problem of launching a vehicle into orbit was simulated 
for flight restricted to the equatorial plane. Difficulties arose at low altitudes 
due to the large and highly nonlinear aerodynamic forces. While the expan- 
sion method gave reasonable results, the use of aerodynamic pulse functions 
in order to reshape the zeroth-order trajectory was vital to matching the opti- 
mal trajectory. Thus it is essential that the zeroth-order path, upon which the 
higher-order corrections are based, resembles the optimal solution such that the 
assumed perturbing effects are indeed small. Based on the difficulties caused 
by the aerodynamics and the modelling of these aerodynamics a few sugges- 
tions are offered. First, the modelling of atmospheric and aerodynamic terms 
should be adequate well beyond the domain of the optimal solution. This is 
especially necessary if in some manner the zeroth-order trajectory significantly 
deviates from the optimal trajectory. It is also suggested that this technique 
would work better with a symmetric version of the ALS vehicle configuration 
by eliminating the irregular behavior of the Hamiltonian. The results of this 
research showed that the idea of using perturbation theory to perform real-time 
on-line guidance is a valid one. The improvement in computational speed and 
effort over the generation of optimal solutions is evident. Still, the technique as 
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presented here was not designed for computational efficiency. To that end, use 
of parallel processing in the integration of the quadratures is proposed. This 
has the potential of decreasing the computational time even further. Lastly, 
one of the goals for deriving an on-line guidance scheme is to provide abort op- 
portunities or instantaneous changes in the terminal destination. It should be 
remembered that this method always provides a nominal path which satisfies 
the terminal constraints while an improvement in the performance is obtained 
by the first-order correction terms. Because of the robustness of the solution 
to the zeroth-order analytic two-point boundary problem, in-flight aborts can 
be easily included. 

More sophisticated modelling, such as including an oblate Earth model 
and wind profiles, can easily be done since these effects can be considered per- 
turbations and included in the problem in the higher-order correction terms. 
The result would be to integrate some additional quadratures. The technique 
can also be extended in a straightforward manner to flight in three-dimensions 
in order to reach a point in space. See appendix [A] for the zeroth-order ana- 
lytic solution. This is done through the addition of the out-of-plane equations 
but with an accompanying increase in the complexity of the problem. It is ex- 
pected that the three-dimensional solution will increase the payload available 
at orbital insertion due to the benefits of the rotational effects of the Earth. 
The inclusion of a dynamic pressure point inequality constraint is also feasible. 
In other ALS studies [4, 12] it has been shown that since the rocket cannot be 
throttled, the vehicle only touches the dynamic pressure constraint at a point. 
It is suggested here that this point inequality constraint can be included in the 
analytic solution of the zeroth-order problem as an interior point constraint 
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[see appendix C]. This new zeroth-order trajectory would avoid the large aero- 
dynamic correction terms found in the unconstrained optimization problem. 
These large correction terms are due to the aerodynamic forces encountered 
when flying through the region of high dynamic pressure. The solution to a 
zeroth-order trajectory including a dynamic pressure constraint represents the 
most important step in the evolution of this method for implementation as a 
real-time, on-line guidance scheme for the launch problem. 

Research is in order, on the inclusion of variable state and control in- 
equality constraints in the expansion of the Hamilton- Jacobi-Bellman equation, 
so as to generalize the class of optimization problems amenable to expansion 
techniques. Future studies are required on the general properties of the validity 
of the asymptotic expansion of the HJB equation. For example, the question 
of whether or not the asymptotic expansion is uniformly convergent remains 
to be answered. This is especially true in light of the use of the ad hoc aerody- 
namic pulse functions. Also, the approximate optimal guidance scheme must be 
made robust with respect to parameter variations and stochastic disturbances 
in order to be implemented as a real-time on-line guidance scheme. Possible 
solution methods [21, 22, 23] have been proposed to handle the more realistic 
situation of a nondeterministic environment. The use of the Hamilton- Jacobi- 
Bellman equation is particular advantageous under these circumstances since 
a stochastic version of the equation exists. With the zeroth-order trajectory 
providing full state information, the best solution in the presence of random 
disturbances should be obtainable. 


Appendix A 


Zeroth-Order Solution for Three-Dimensional Flight 


The analytic zeroth-order solution is derived once again by a transfor- 
mation of coordinate system. A canonical transformation from the wind axis 
to the rectangular or local horizon coordinate frame allows the zeroth-order 
problem to be solved analytically. The solution is in closed form up to some 
constants that can be determined numerically. By making the transformation 


u — V cos 7 cos x 

(A.l) 

v = V cos 7 sin x 

(A.2) 

w — — V sin 7 

(A. 3) 

the zeroth-order equations of motion in a cartesian coordinate frame become 

X = u 

(A.4) 

Y — v 

(A. 5) 

h = —w 

(A. 6 ) 

T 

ii = — cos 0 V cos ip 

(A- 7) 

m 

v = — cos^„sin^ 

(A- 8 ) 

m 

w = sin 0 P + g 3 

m 

(A.9) 


m = —crT 
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where the Thrust pitch attitude, O v = a 4- 7, and the Thrust yaw angle, xp = 
X + become the control variables for the zeroth-order solution. 

The optimization problem to zeroth-order is solved by the Hamilto- 

nian 


H = \xu + X yv — A uw + A u — cos 0 P cos xp + 


m 


r r 

X v — cos 0 P sin xp + X w ( sin 0 V + g s ) 

m m 


(A. 10) 


The zeroth-order control laws determined by the optimality conditons are: 


Thus, 


T 

Hip — — cos 0 P { X v cos xp — A u sin xp) = 0 

T T 

Ho„ = -(A u cos xp + X v sin xp) — sin 0 V A,„ cos 0 V = 0 

m m 


(A. 11) 
(A. 12) 


taxixp = 
cos xp = 

sin xp = 


K 

A„ 


A u 


\Au + 

Au 


•JK + ^ 

and using these forms for the cosine and sine of xp, 0 P is 

A,., 


tan 0 P = 


cos 0 P = 


sin 0 P = 


\Au + ^u 


yj^t + Ag 

\Au + 


(A. 13) 


(A. 14) 


Propagation of the Lagrange multipliers by X x = -//J gives 

Ax = 0 

Ay = () 

A/, = 0 
A u = — Ax 

An = — Ay 

A w — A h 

with boundary conditions given by X f - C xJ 
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(A. 15) 


A x(Tf) = i/ x , Ay(r/) = u Y , \ h ( Tf ) = u h , A u ( Tf ) = u u , 

a v( T f) = i>v, A w (r f ) = i/,„ 

where v x , vy , v h , u u , i> v , i/ w arc unknown I.agrange multipliers associated with 
the terminal constraints. The solutions to the adjoint differential equations are 

Ax = vx 

Ay = l/y 

A/> = i>h, 

K = C u - Ax(r - r 0 ) 

An = C v — Xy{t — 7"o) 

Aw = C w + X h (r - T 0 ) (A. 16) 

The equations of motion can be integrated by changing the independent vari- 
able from time to mass and using the mass equation to substitute mass for r. 
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Therefore, the Lagrange multipliers become 


Au 

= C u + X x ^ 

al 

(A. 17) 

A„ 

= C v + \y ^ 

(A. 18) 

X w 

7^ X m 

'U> ^r]-, 

(A. 19) 

X 2 u + Xl + Ai 

= cm 2 4* bm + a 

(A. 20) 


where 


c = 

Ax + Ay- -f X 2 k 

{oTY 

(A. 21) 

b = 

4p(X xC u + X Y C V - X h C w ) 
al 

(A. 22) 

a = 

cl + cl + cl 

(A. 23) 

= 

c u - Xx — 

(A. 24) 

= 

r \ m ° 
a-A v — 

(A. 25) 

= 

c w + x h ^ 

(A. 26) 


and the state equations become 


du 

X u 


(A-27) 

dm 

amy/cm 2 + bm + a 


dv 

x v 


(A. 28) 

dm 

amy/cm 2 + bm + a 


dw 

Au, 

9s 

(A. 29) 

dm 

amy/cm 2 + bm + a 

oT 

dX 

u 


(A. 30) 

dm 

~oT 


dY 

V 


(A.31) 

dm 

~oT 


dh 

w 


(A. 32) 

dm 

of 
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Note that c > 0, a > 0 and A = 4ac — b 2 > 0 since 

A = - ■ — t [(AxCu — AyC u ) 2 + (A xC w 4- A/,C U ) 2 + (\yC w + A/,C„) 2 ] (A. 33) 

{al y 

Therefore, the integrations can be obtained from standard integrals. After some 
manipulation, solutions for the states are obtained in terms of the unknown 
constants as 
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4-rn 


Ay- , ,b 

^Piy7c (m + ^ 

c„ 


sinh-' fet* ) - slnh -' 


_1 / 2 cttto 4 - fa' 


+ 


a(aT)y/a 
C v 


sinh 


V ^ 

-i / 2 a + fam \ 


m 


VE 


sinh' 


V 7a 

2a + faroo 


m 0 


v/A 


h = h o — 


<j{<jT)y/c 

Ay 

cr(aT) 2 c 

9s 


. l {2cm + b\ . _j [2cmo + fa' 
sinh j= — — sinh 7 = — 

l j l *5 , 


yi ct 7 Zq + fam 0 4 - a — \/ cm 2 4 - fam 4 - a 


(A. 38) 


2(<rT) 2 

A a 


(m — m 0 ) 2 4 - 


(m - mo) 
aT" 


-w 0 


a(^VS (m+ ^ ) 


. / 2cm 4- fa' 

sinh p= — 

V vS , 


— sinh' 


/ 2 cmo + fa' 

v~ 7 a - , 


— m 


a(<jT)y/a 
C w 


. . / 2 a 4- fara\ . , /2a 4 fam 0 

sinh ■=- — sinh =- 

\ mv A / \ rrioV A 


o(oT)yJc 

A/i 

cr(<jT) 2 c 


sinh- 1 ^2cm + 6^ _ stah _. 


l V5 


2 cm 0 4 - fa' 
, \/A , 


y/ crrifl 4 - fam 0 4 - a — 'Jerri 2 4 - fam 4 - a 


(A. 39) 


There are seven unknown constants that are to be determined 
m/, C u , Cu, C™, Ax, Av, A*. These unknown constants can be found using the 
initial and final states which arc known, the six state equations above, and the 
transversality condition for the Hamiltonian. 


These equations are valid for arcs before and after staging occurs. To 
determine a point on the trajectory after staging, given initial conditions before 
staging, the Wei ers trass- Erdmann corner conditions can be used to relate the 
Lagrange multipliers before and after staging and thus link the two subarcs. 
Since the states are assumed continuous across the stage time and the change 
in mass is a known fixed quantity, all the Lagrange multipliers are continuous 
in time. The Return Function is thus continuous and constant across the stage 
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time. But since the analytic state equations arc derived using mass as the 
independent variable, the equations for the states and Lagrange multipliers 
change across the stage. The form of the equations is the same but the initial 
values for the states in the equations are now replaced by the values of the 
states at staging before the discontinuity in mass. The mass flow rate, crT, will 
change and the initial time associated with the Lagrange multipliers becomes 


the stage time. Therefore, 


A u (£) — ( A u (L*t) — Ax 


m st 2 


(tT s 


12 


— C u T A; 


rn 


a T Si 


a 


A v (i) = A v {t sl ) — \y 


m st 2 


<tT< 


si 2 


— C v + Ay 


m 




sV2 


A w (0 — ( A^(i s t) + Ah — ; 

o t SL2 




A* 


m 


&T sl 2 


+ Ax 


m 


aT< 


st 2 


-4- A y 


m 


oT 5 


st2 


— A k 


m 


0Ty 


st2 


Thus the constants a,b,and c in the analytic solution have the same form but 
the term C after staging becomes related to C before staging. 


v. 


C u — Ax 


— A y 

C w + A/ t 


' m st2 

m st i 

o r l\ a 

crT sa 

m. s t2 

m 3l i 

oT sl 2 

crT sl i 

m st2 

m sl , 

oT sl2 

vT 3t i 


) 

) 

) 


Through a coordinate transformation back into the wind axis the an- 
gle of attack can be determined. Higher-order terms in the Lagrange multipliers 
can be found from expansion of the dynamic programming equation. The Opti- 
mal Return Function can be determined by integrating the perturbation terms 
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along the zeroth-order trajectory while taking into account variations in the 
final time due to changes in the initial states. Taking the partial derivative of 
the Return Function with respect to the initial states determines the Lagrange 
multipliers. Expanding to first-order the control and Lagrange multipliers in 
the control law determines the first-order approximation for the controls. 

A.l Zeroth-Order Coordinate Transformation 

The analytic solution for the zeroth-order problem has been found 
in the Cartesian coordinate system but the equations of motion of the full 
system which includes the aerodynamic forces are written in the wind axes 
system. To derive the zeroth-order control and the first-order correction to the 
controls the transformation of coordinates and especially the transformation of 
the Lagrange multipliers must be known. The rotation from the wind axes to 
the local horizon frame is done by a canonical transformation. 

u = V cos 7 cos x 
v — V cos 7 sin x 

VJ = -V'sin7 (A. 40 ) 

A necessary and sufficient condition for the transformation to be canonical is 
that the Hamiltonians be equivalent. 

H lh — A xd.X + XydY -f- A hdh 4- X u du + X V dv + A W dw 
Hwind = A gdO + X ^dif) -I- A h dh + XydV + X^d^ + A X dx 
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For the two reference frames, h and m are the same. Also, for the model 
X = r c 0 and Y = r e (p. Thus, 

X 0 = r e Xx, A* = r e Xy , and X h is the same 

In order to equate the Hamiltonians, the Jacobian of the transforma- 
tion from the wind axes to the local horizon axes is computed. 

X u du + X v dv + X w dw = XvdV + A 7 d7 + A x d\ 


So, 


Av 


A. 

= 

A x . 



Du 

Ov 

Dw 

av 

DV 

DV 

Du 

Dv 

Dw 

Dy 

D'y 

Dy 

Du 

Dv 

Dw 

Dx 

Dx 

dX 



c 

1 


A v 


1 

£ 

i 


and thus 


Av 
A 7 

— 

1 



cos 7 cos x cos 7 sin x -sin 7 
-V sin 7 cos x -Fsin7sinx -V cos^ 
- V cos 7 sin x V cos 7 cos X 0 


A u 

A v 

A„, 


Therefore, 


Av 

A 7 

A x 

A 0 

A<^ 


A u cos 7 cos x + A„cos7sinx - Au,sin7 (A. 41 ) 

— F(A U sin 7 cos x + A„sin7sinx + A.VCOS7) (A. 42 ) 
\/cos7(— A u sin x + A v cosx) (A. 43 ) 

r e Ax 
r e \y 


and 


V 

tan x 


\/u 2 + t / 2 + tu 2 

V 

u 


w 

V 


sin 7 













Appendix B 


Canonical Transformations 


The use of the canonical transformation of section 4.2 is essential in 
finding the closed-form solution to the zeroth-order problem because canoni- 
cal transformations preserve the Hamiltonian form of the equations of motion 
in the new set of variables. A more thorough discussion of canonical trans- 
formations than what is presented here is in [24]. To transform between the 
generalized coordinates and generalized momenta or Lagrange multipliers ( q , p) 
of one system to new variables (Q, P) of a new system, a set of transformation 
equations linking the two systems must be known. This link between the two 
systems can be derived from the generating function S(q,Q,t). 

4 S(q,Q,t) = L(q,q,t) ~ L(Q,Q,t) (B-l) 

(it 

where L I — ■ V is the Lagrangian of the respective system. Let the transfor- 
mation equations be of the form 

Q = Q(q,p,t) P=P(<l,P,t) (B.2) 

with Hamiltonians associated with each set of variables such that the Hamil- 
tonian equations are satisfied, i.c. , 

H(q,p,t) = - £(g,<M) 

(PL 

dq 


Pi 
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(B.3) 

(B.4) 
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(13-5) 

(B.6) 

(B.7) 


<h = 

K{Q,P,t) = 

P. = 


Qi 


dir 

dp i 

J2 PxQx - L(Q,Q,t) 

dL 

0Q 

OK 

m 


(B.8) 


By solving Eqs. (B.3)-(B.6) for the Lagrangians and substituting into Eq. (B.l), 
the difference between two differential forms can be obtained as an exact dif- 
ferential. 

^ ZP'd-Qi ~ II dt - YL PidQi + Kdt = dS (B.9) 

This is the sufficient condition for a canonical transformation between the old 
variables (q, p) and the new variables (Q, P). 

The generating function can be written as a total differential of the 

form 

, c ^ dS , dS ^ dS J 

iS = ^^' + ^0Q^ Q ' + K dt ‘ B10 > 

Equating like terms of the differential dS yields 


P, = 


dS_ 

dq % ’ 



K = II + 


dS 

dt 


(13-11) 


A simplified form of these equations can be obtained. Assume that 
time is not changed in the transformation from one system to the other system. 
Therefore, t is an independent parameter and the value of dt is set to zero. Also, 
define a new function equal in value to the generating function but expressed 
in a different form 


tffa.P.O = S(q,Q,t) 


(13-12) 
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Therefore, the variation ol Hq. (B.l) is rewritten ;is 

Z pM - Z ' w, = ^ ( B - 13 ) 


The equations 


dip 

% 

dip 

dpi 


a-Z^ 


5^ 


- _ v P 

r J fa 


(B.14) 
(B. 15) 


are obtained by expressing the variations tiiy and P^SQ, in terms of the old 
variables. The new Hamiltonian is determined by integrating the expressions 


presented above to obtain ip, and then calculating 


K = II + 


dip 

~di 


+ 


y P .?8i 

^ 1 dt 


(B. 16) 


For a homogeneous canonical transformation the generating function 


S or ip is identically zero and thus dip is an exact differential and equal to zero. 
For a point canonical transformation as used in section 4.2, the transformation 
equations Q t = Q l (q) are functions only of the generalized coordinates of the 
old system. The functions Q t (q) represent a full set of independent functions, 
therefore, the old variables q can be expressed in terms of the new variables Q 
using these same functions. This implies that the Jacobian determinant is not 
zero, 


■ dQj . _ d{Qu Qjh •••> Qn) _j_ q (B.17) 

dqj d{q\,q 2 ,--,qn) 

Since the transformation equations (Q, ip) do not contain time, the Hamiltoni- 
ans of the two systems are equal (K = //). This result is obtained by the use 
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of Eq. (B.16). The relationships expressed in Eqs. (B. 14)-(B. 15) become 


0 

0 


p,-Zr,i Q ’ 




j % 

OQj 


dpi 


(B.18) 
(B. 19) 


The Jacobian determinant of the new variables Q in terms of the old set of La- 
grange multipliers p must be zero for a set of P’s not all identically zero and is 
also required since the point transformation Q = Q(q) was independent of the 
Lagrange multipliers p. The matrix ^ is used to determine the relationship 
(Eq. (B.18)) between the Lagrange multipliers of the two systems. The trans- 
formation between the wind axis coordinates and the Cartesian coordinates of 
section 4.2 is a canonical transformation as can be verified by the use of the 
canonical transformation equation (Eq. (B.13)). 



Appendix C 


Point Inequality Constraints 


The inclusion of the a point inequality constraint on the dynamic 
pressure is discussed in this section. Due to structural load limits imposed 
on the ALS vehicle, the optimization problem must include a dynamic pres- 
sure inequality constraint. For the unconstrained optimization problem that 
is presented in this research, the correction terms to the zeroth-ordcr solution 
become too large near the region of maximum dynamic pressure. For a rocket 
incapable of throttling, the optimal trajectory will not include a subarc on the 
boundary of the dynamic pressure constraint but instead will only touch the 
constraint at a point [4, 12]. This result would seem to indicate that the dy- 
namic pressure inequality constraint can be handled in the same manner as 
an interior point constraint. Therefore, the Ilamilton-Jacobi-Bellman equation 
can be split into subarcs before and after satisfying the interior point constraint. 
Let the optimal return function be 

OO 

P{x,L) = '*Tc'P l (x, t) (C.l) 

t=() 

The point interior equality constraint is 

N(y(tx)) = 0 (C.2) 

where the constraint is a function of the states y at the time ii and the states 
are assumed continuous at the constraint. The system dynamics are defined 
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by 


y = f{y,u,t) (c.3) 

and therefore the Hamiltonian is II — P x f ■ Across the equality constraint 


Px(]J,t\ ) — Px{y, t\+) + NN x (y,t i) (C.4) 

Px(y,ti-)f(y,u,ti-) = P x (y,t ] +)f(y,u,t l +) (C.5) 


where El is the Lagrange multiplier used to adjoin the constraint to the per- 
formance index. These equations are the corner conditions derived using the 
calculus of variations [9]. From the solution of the Hamilton-Jacobi-Bcllman 
first-order partial differential equation, the higher-order terms of the expansion 
of the optimal return function are 


Pi(x,t) = -J^ Ilidr — J^+Ilidr i= 1 , 2 ,... (C.6) 

Recall that the integration is performed along the zeroth-order trajectory. The 
partial of the return function with respect to the initial state x is 




V + ‘^ yul,+) ^ - l' + Tit dT (C7 > 

This equation determines the correction terms to the Lagrange multipli 


dx 
dt i 


lers. 


Notice that the variation in the time at which the equality constraint is satisfied 
is explicitly taken into account in the correction terms. Substituting the partial 
of the expansion of the return function (Rq. (C.l)) into the corner condition of 


Eq. (C.4) produces 


£ *i-) = £ JPiAv, ti+) + nN x ( y , u) 


t=0 


t=0 


(C.8) 
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The solution for the Lagrange multipliers is then determined by collecting the 
coefficients of like powers ol the expansion parameter, c. Iherefore, for the 
zeroth-order term of 2 = 0, 

Po x {yA\ — ) — ^o x (?y> 1 1 +) + UN x (y> t\) (C.9) 

and for higher-order correction terms 

PiAy>tx-) = PiAv>t\+) 1 ^ 1 ( c - 10 ) 

This result implies that the higher-order correction terms for the Lagrange mul- 
tipliers are continuous across the equality constraint. The jump discontinuity 
in the Lagrange multipliers is completely taken into account by the zeroth- 
order term. The higher-order terms of the expansio' of the return function 
are thereby continuous across the corner. I he continuity of the Hamiltonian is 
ensured by substitution of the partial of Eq. (C.l) into Eq. (C.5), which results 
in 

oo 00 

yV/U^i-)/(ZMO>-) = £e t a r (2/,ii + )/(2/,Mi+) (0.11) 

t=0 >=° 

Thus, by collecting terms in like powers of c, the equation 

Pi x (x,t\-)f(y,u,L\-) = Pi x {x,t\+)f{y,u,ti+) (C.12) 

is obtained. This condition is just the continuity of the expansion terms of the 
Hamiltonian, i.e. , //,(?/, u, P x , 1 1 — ) = P x ,t\+). 

The dynamics of the system are the same after meeting the point 
constraint as they were before; the constraint. Ihus, the analytic solution of 
the states as derived previously is still valid but with a change in the La- 
grange multipliers at the dynamic pressure constraint. The relationship be- 
tween the Lagrange multipliers across the constraint given by Eq. (C.9) can be 
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used to link the subarc of the zeroth-order trajectory that occurs before the 
constraint equality is met to the subarc that occurs after the constraint equal- 
ity is met. These equations are additional conditions that are used to solve the 
constants associated with the two-point boundary value problem. Since the 
zeroth-order problem totally accounts for the dynamic pressure equality con- 
straint N(y,t\) = 0, the first-order correction term to the zeroth-order term 
should be small and the asymptotic expansion should remain valid near the con- 
straint. Therefore according to the method of characteristics, the zeroth-order 
trajectory can be used as the characteristic curve to determine the higher-order 
correction terms. In contrast to the unconstrained optimization problem, this 
zeroth-order trajectory should be close enough to the optimal solution such 
that only small corrections to the zeroth-order guidance law are necessary. 


Appendix D 


Analytic Partial Derivatives for Zeroth-Order Solution 


The partial derivatives with respect to the general initial state x de- 
rived in Chapter 5 for the first-order correction terms are presented here in 
their explicit form for the initial states, Vo and 7 o- fl ic equations derived are 
for the second stage subarc. The initial velocity components expressed in terms 
of the wind axis states are 


iiq — Ci cos 7o, u.'o — Vo sin 70 


(D.l) 


and therefore the partial derivatives with respect to the initial velocity and 
flight path angle are 


du 0 
OV 0 
du 0 

<? 7 o 
dw 0 
dV Q 
du: 0 

d 70 


cos 70 


-Vosin7o 


- sin 70 

- V 0 cos 70 


(D. 2 ) 


The partials of the analytic zeroth-order state equations are then expressed as 


du 

Wo 


cos 70 


C u 


r y/a 


<7 y/oL 


'8C U C u da ' 
K dV 0 2a < 9 Vo 

d% 1 3^2 (mo) 


sinh 1 — sinh 1 ^(^o)] 


7 l 4 - 31 (m) ,,v « v ,r l + !lV '‘ 


(D. 3 ) 
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du 

(ho 


— — Vosin7o 7= [sinh 1 $ 2 (m) — sinh 1 S 2 


a 


CTy/a L 
1 



dC u 

(ho 


dQo 


d$s 2 (m o) 


L\/l+^(m) ho + ^ (mo) d l0 


C u da \ 

2 a d 7o J 

(D.4) 


<9w 

dV 0 


dw 
d~i o 


sin 70 - 
C w 


1 


<7v/a 

1 


sinh 1 ^s 2 (m) — sinh 1 

dS 2 1 



dC w C w da 


d V 0 2 a dV 0 


. d%(m 0 ) 

‘'v'a |y l + 3*( m ) <9Vo + 3 |( mo ) SV 0 


1 


a 2 T y/c 


A* 



■ V Q cos 70 7= [sinh -1 3 2 (m) - sinh -1 S 

ay/a L 



dC w 

d'Yo 


' W 


ay/a 

1 


(D-5) 

C w da \ 
2 a (97 o / 


1 (h >2 _ 1 d^s 2 (m 0 ) 

\/l + ^(ra) y^l + S 2 (m 0 ) *^7o 


2 Ty/c 

Aa 


(7 


2T 



(D.6) 


<9/l 


- sin 7 q 

m 


aT y/a 

C w 


(m - mo) 

aT 

sinh- 1 0 2 (m) - sinh- 1 3 2 (m 0 )l (%£ - £=J£A 

J \ dvQ 2 a dVo J 


-m- 


1. 


00, 


d%(m o) 


aT yfi[ ^J^T%Um) ^v o yf\ + ^( mo) dVo 
^ m 2a{aT) 2 CT 2 [ s ' n ^ ^i( m ) ~ sinh — 
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dh 

d'yo 


|sinh 1 ^i(rn) - sinh 1 Gi(ra 0 )] 


'em? 4- bm 4- 


o{pT') 2 \fc 

1 t 

a{aT) 2 c [ 


2a(aT) 2 c ^ \J cm? 4- bm + a 

(m - mo) 

■ sin 7 o = — 


— T~Z 1 ( dXh A* dc 

cm& + bmo + a\\—-- 


dV 0 c dV 0 


+ gfe 


__2 dc i db | da 

m 0 3 Vq m ° 9V 0 dVp 

Jcml + bmo 4- a 


771 r • I -1 rx. / \ • l 1 <x / Cm da 

■rfvsi smh 32(m) - s '" h 

C w [ 1 d3 2 1 dQ 2 (mo) 

— yrfl — — — ■" — - 

vTs/a y/\ 4. 'Csl(m) ^7o y l + 0^(m 0 ) ^7° 

I 5 '" 1 ’"' a ' <m) " Sintr ' °' (m0) l Wo 

~WWTc I sinh "' a,(m) “ sinh "' a,(mo) ] W 


<y{aT) 2 c 


[sinh 1 li| (rn) — sinh 1 ^i(mo)] 

, — / — : ] ( d\ h A h dc 

Jem? 4- bm 4- a — J cttlq 4- bm 0 4- a I — — — 


L.UL UUC ~T~ U. Y K ' ,,L 0 ' w,IL i) J ^ Q 

2 9c I „ <>b , <>o 2 dc . db , da 1 

771 <>70 + m O70 + <>7n m O»ro + 7770 <>ro + <>7o /r 


2 a(aT) 2 c \J cm? 4- bm 4- 


>70 ^ ^<>70 ^ <>70 

crrio 4- brrio + a 


The partial derivatives of the eonstants a,b,c, and C w used to express 
the analytic state equations arc 


da 

Wo 

da 

d'yo 

db 

dV 0 

db 

d'yo 


on < ^ u , 9/^ dCm 

u dv 0 w dv 0 
2 C^ + 2 cJ^ 

d'yo o 7o 


(D.9) 
(0. 10) 
(D. 11) 
(D.12) 
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5c 

dV 0 

dC w 

dV 0 

dc 

dlo 

dC w 

dv 0 

8C W 

d7o 


2 d\ h 

(crT) 2 5 Vo 
dC w mo 5 A* 
5V& + <rT dV 0 
2 5A„ 

{ary 5 7o 
dC w m o d\ h 
~dv { o’ + 7r5v^ 
dCw rno d\ h 
dj 0 + crT c>7o 


Also, the derivatives of the discriminant function A = 4ac — 6 2 are 


d A 
5K 0 
5A 

$70 


4a 


5c 


da 

+ 4c—— - 26 


db 


dV 0 d V 0 dV 0 

dc da db 

4a- I- 4c— 26— — 

57o 57 q 57o 


(D. 13) 

(D. 14) 
(D.15) 

(D.16) 

(D.17) 

(D. 18) 
(D. 19) 


The arguments of the inverse hyperbolic sine function are defined in 


Eq. (5.19) as 


3i (m) 


2 cm 4- 6 

^7^ 


9 2 (m) = 


2a 4 6m 

mV A 


The resulting partial derivatives of the arguments are 


53j 

5V^ 


59 1 

57o 


59 2 

0V O 


1 

71 


i 

7a 


dc ( n , 69 i n 56 0 Sq 5a 

( vs’w. ( VE } dV 0 c vSav 0 


+2c 


5m 

5K) 


(D.20) 


m 


7a 


a9i .5c t 69 1 56 9 X 5a 

+2c^ 

5 70 

c9 2 5a 69 1 . 56 

2( " m 73 iw> +m 1 + 71 W 


( 0 . 21 ) 


9 2 5c 2 dm 

~ ma VKdV 0 ~ m a Wo 


(D.22) 
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dS 2 

dlo 


_1 

myfK 


. cQ 2 s da . 

2(1 - m—^=)— 4- m( 1 + 


VA d'yo 


b% db 
y/K d'yo 


—ma 


5 2 dc 
y/Kd'io 


2 dm 
m 1 d 7o 


(D.23) 

(D.24) 


Note that unless the partial deriviatives are evaluated at the terminal manifold 
the partials and are zero. Using the partial derivative chain rule for a 
trignometric function, the partials of the inverse hyperbolic sine function are 
obtained. 


w ( sinh '' 3i ) 
Wo ^ 


1 ash 

yrTs? dv Q 

i ds 2 

yr + q 2 a7 ° 


(D.25) 

(D.26) 


Therefore, all the partials derivatives needed to evaluate the partial derivatives 
of the analytic state equations along the zeroth-order trajectory depend on the 
eight constant partial derivatives These 

partial derivatives are functions of the solution to the two point-boundary value 
problem. Therefore, they are constant when integrating the forcing function R\ 
from the initial to the final conditions but they change as new initial conditions 


are given when the guidance scheme is used in feedback. 
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